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Abstract 

Estimation  theory  is  applied  to  a  physical  model  of  incoherent  polarized  light 
to  address  problems  in  polarimetric  image  registration,  restoration,  and  analysis 
for  electro-optical  imaging  systems.  In  the  image  registration  case,  the  Cramer- 
Rao  lower  bound  on  unbiased  joint  estimates  of  the  registration  parameters  and  the 
underlying  scene  is  derived,  simplified  using  matrix  methods,  and  used  to  explain  the 
behavior  of  multi-channel  linear  polarimetric  imagers.  This  discussion  is  expanded 
to  include  biased  estimators  up  to  the  point  where  the  results  become  registration 
algorithm  specific. 

In  the  image  restoration  case,  a  polarimetric  maximum  likelihood  blind  decon¬ 
volution  algorithm  is  derived  and  tested  using  laboratory  and  simulated  imagery. 
The  angle  of  polarization,  polarized  and  unpolarized  components  of  the  scene,  and 
channel  point  spread  functions  are  jointly  estimated  using  this  approach.  The  results 
of  this  estimation  are  combined  to  unambiguously  determine  the  scene  state  of  linear 
polarization. 

Finally,  a  principal  components  analysis  is  derived  for  polarization  imaging 
systems.  This  analysis  expands  upon  existing  research  by  including  an  allowance  for 
partially  polarized  and  unpolarized  light.  From  these  results,  a  means  of  detecting 
polarizing  objects  under  weakly  polarizing  obscurations  is  proposed  by  discarding 
the  principal  component  channels  that  are  insensitive  to  temporal  fluctuations  in 
scene  polarization. 
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I.  Introduction 

This  dissertation  describes  research  in  three  areas  of  passive  polarimetric  im¬ 
age  processing:  image  registration,  restoration,  and  optimized  analysis.  This 
chapter  provides  an  overview  of  polarimetric  imaging  in  remote  sensing,  describes 
the  problems  addressed  in  this  research,  and  provides  an  overview  of  this  document 
as  a  whole. 

1.1  Polarimetric  Imaging 

Polarimetric  Imaging  (PI)  is  a  form  of  remote  sensing  in  which  an  object  or 
event  is  characterized  in  terms  of  the  polarization  state  of  its  reflected  radiation. 
Polarization,  which  is  a  fundamental  property  of  electro-magnetic  radiation,  will  be 
defined  in  detail  in  Chapter  II.  The  purpose  of  this  section  is  to  introduce  and 
motivate  polarimetric  imaging  as  a  remote  sensing  discipline. 

Though  less  common  than  panchromatic,  multispectral,  and  radar  imaging, 
polarimetry  has  been  applied  successfully  to  a  number  of  remote  sensing  problems. 
Astronomers  have  used  PI  to  characterize  the  surface  of  planets  and  their  atmo¬ 
sphere,  to  determine  the  surface  properties  of  comets  and  asteroids,  and  to  investi¬ 
gate  novae  and  the  interstellar  medium  [4,20].  In  the  earth’s  atmosphere,  attempts 
have  been  made  to  characterize  clouds  and  atmospheric  aerosols.  Additionally,  PI 
has  been  applied  to  earth  surface  characterization  problems  such  as  hydrology  and 
oceanography  as  well  as  inventory  and  species  identification  problems  in  forestry 
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and  agriculture  [8,47].  The  reader  may  also  be  interested  in  applications  of  PI  in 
machine  vision  [27]  and  imaging  through  scattering  media  [44],  Finally,  additional 
applications  are  discussed  in  texts  dedicated  to  the  subject  of  polarization  [5,42]  as 
well  as  general  optical  texts  [3,16,31]. 

The  military  applications  of  polarimetric  imagery  include  targeting,  reconnais¬ 
sance,  and  intelligence.  The  targeting  and  reconnaissance  applications  are  straight¬ 
forward:  man-made  materials  tend  to  produce  more  highly  polarized  reflections  than 
their  surroundings;  ergo,  polarimetry  is  a  queueing  aid.  In  addition,  adversary  cam¬ 
ouflage,  concealment,  and  deception  techniques  are  becoming  increasingly  sophisti¬ 
cated  and  include  those  with  spectral  and  radar  defeating  properties;  PI  is  an  addi¬ 
tional  discriminator  that  can  be  used  to  defeat  these  techniques  separate  from  or  in 
addition  to  other  sensing  modalities.  PI  has  also  been  shown  to  be  a  capable  detector 
of  smoke  and  rocket  plumes  [7].  Besides  intrinsic  value  as  a  queueing  technique,  PI 
can  potentially  be  used  as  a  material  identification  and  surface  characterization  tool 
especially  when  used  in  conjunction  with  spectral  analysis.  As  alluded  to  previously, 
PI  can  be  used  to  characterize  asteroids  and  it  has  been  suggested  that  PI  can  be 
used  to  locate  and  characterize  space  debris  [47].  Finally,  Strong  makes  an  argument 
for  PI  as  an  enabler  for  space  situational  awareness  [43]. 

The  literature  cited  in  this  section  is  provided  as  an  overview  of  the  existing 
work.  To  facilitate  understanding,  a  more  complete  literature  review  on  each  of 
the  addressed  aspects  of  polarimetric  imaging  is  included  in  the  beginning  of  each 
chapter. 

1.2  Research  Contributions 

This  work  is  an  estimation  theory  approach  to  processing  polarization  imagery. 
The  references  in  the  preceding  section  solve  problems  using  polarization  imagery 
and  from  a  physics-based,  phenomenological  perspective.  In  other  cases,  such  as 
the  works  listed  in  the  introduction  to  Chapter  III,  estimation  problems  (e.g.  image 
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registration)  are  treated  as  equivalent  to  their  polarization  insensitive  counterparts. 
(Noteworthy  exceptions  to  this  rule  include  [43,51].  Consequently,  estimation  re¬ 
search  into  polarization  imagery  is  ripe  with  possibilities,  even  in  cases  that  have 
been  exhausted  in  traditional  image  processing  literature. 

Three  polarimetric  imaging  processing  problems  are  addressed  in  this  disser¬ 
tation:  registration,  restoration,  and  optimized  analysis.  Polarimetric  imaging  is 
inherently  a  multi-channel  discipline.  A  channel,  in  this  context,  refers  to  an  imag¬ 
ing  path  that  is  most  sensitive  to  a  specific  state  of  scene  polarization.  Reliable 
channel-to-channel  registration  is  therefore  a  prerequisite  for  exploitation  of  polari¬ 
metric  data  and  worthy  of  further  study.  Image  restoration,  which  is  the  process  of 
estimating  the  true  target  scene  in  the  presence  of  a  corrupting  transmission  media 
and  noise,  is  often  studied  in  traditional  intensity  imaging  to  improve  sensor  perfor¬ 
mance  without  increasing  cost,  size,  and  complexity;  here,  these  concepts  are  applied 
to  polarization  imaging  to  achieve  the  same  advantages.  And,  at  the  end  of  the  pro¬ 
cessing  chain,  the  multichannel  nature  of  these  data  place  an  unnecessary  burden  on 
the  analyst  in  cases  where  some  of  the  channels  contain  no  useful  information.  If  the 
unnecessary  information  can  be  reduced  by  eliminating  channels  (or  combinations  of 
channels)  using  prior  knowledge  then  a  tangible  savings  to  the  analyst  is  achieved. 
An  attempt  to  address  each  of  these  points  is  laid  out  in  the  following  three  chapters. 

In  Chapter  III,  the  effects  of  polarization  diversity,  channel  noise,  and  regis¬ 
tration  errors  on  polarization  image  estimation  are  considered  in  the  context  of  the 
Cramer- Rao  Lower  Bound  (CRLB).  The  bound  is  derived  for  the  A-charmel  case  in 
the  presence  of  random  translational  registration  errors.  The  bound  is  then  used  to 
study  three  and  four  channel  polarization  imagers  and  the  special  case  of  N  misregis- 
tered  polarization  insensitive  images.  The  bulk  of  this  work  was  originally  published 
in  [24],  Some  new  material  on  the  bound,  including  a  discussion  of  external  param¬ 
eter  measurement  versus  joint  estimation  and  a  correlation  based  interpretation  of 
the  CRLB,  can  be  found  in  Appendix  B. 
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In  Chapter  IV,  a  multichannel  blind  deconvolution  algorithm  is  developed  for 
polarization  imagery.  The  algorithm  is  an  instance  of  maximum  likelihood  estimation 
derived  from  expectation  maximization  theory.  This  derivation  and  some  results 
using  laboratory  data  are  published  in  [25].  To  this  work,  the  later  sections  of 
the  chapter  are  devoted  to  new  material  addressing  the  attributes  of  the  estimator 
extracted  via  simulation.  This  new  work  includes  a  comparison  of  the  multichannel 
estimator  to  a  single  channel  estimator  and  a  discussion  of  estimator  bias  with  respect 
to  polarization  angle. 

In  Chapter  V,  a  method  is  proposed  for  optimizing  combinations  of  polarimet- 
ric  imaging  channels.  Optimization  is  achieved  through  identification  of  combina¬ 
tions  of  channels  that  yield  the  greatest  sensitivity  to  polarization  effects,  especially 
for  the  case  of  obscured  targets.  The  remaining  data  are  discarded  as  superfluous. 
These  results  are  obtained  using  principal  components  analysis  (PCA)  and  threshold¬ 
ing  based  on  assumed  temporal  fluctuations  in  polarization  state.  Theoretical  work 
is  described  for  three  and  four  channel  systems  and  simulation  is  used  to  demonstrate 
utility  in  the  four  channel  case.  This  work  was  originally  published  in  [26]. 

Additionally,  several  practical  laboratory  polarimetric  imagers  were  constructed 
and  employed  during  the  course  of  this  research.  Though  not  a  novel  research  con¬ 
tribution  per  se,  this  design  work  is  anticipated  to  be  of  value  to  future  students  and 
researchers.  As  such,  Appendix  D  includes  the  salient  design  details  and  a  discussion 
of  sampling. 

1.3  Scope 

PI  has  applications  in  both  passive  and  active  sensor  systems.  These  systems 
include  those  with  responsiveness  to  emitted  or  reflected  radiation  in  all  spectral 
regimes.  The  scope  of  this  research,  however,  is  limited  to  passive  electro-optical 
remote  sensing  of  incoherent  radiation.  To  avoid  repetition,  the  reader  should  infer 
this  restricted  definition  of  PI  throughout  the  remainder  of  this  document. 
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Polarization  imagers  may  be  comprised  of  one  or  more  channels.  A  single 
polarization  channel,  either  stationary  or  rotating,  can  serve  as  an  aid  for  target 
queueing  [2],  The  two  channel  case,  sometimes  referred  to  as  polarization  difference 
imaging,  is  used  for  target  queueing,  coarse  material  identification,  and  imaging 
through  optically  scattering  media  [41,44],  As  shown  in  Chapter  II,  a  minimum 
of  three  channels  is  required  to  unambiguously  determine  the  state  of  linear  polar¬ 
ization  and  four  or  more  channels  can  sense  both  linear  and  elliptical  polarization 
states.  Since  linear  polarization  is  more  common  than  elliptical  in  remote  sensing  ap¬ 
plications,  the  scope  of  this  research  is  generally  restricted  to  the  linear  polarization, 
three  or  more  channel  case. 

Polarization  imagery  is  extremely  sensitive  to  the  relative  orientations  of  the 
source,  target,  and  receiver.  Assuming  there  is  motion  in  the  target  or  receiver,  the 
validity  of  the  measurements  in  a  multi-channel  polarimetric  imager  have  a  tem¬ 
poral  component.  Throughout  this  research,  all  measurements  are  assumed  to  be 
essentially  simultaneous  on  the  time  scale  of  the  intended  target  scene. 

1.4  Organization 

The  remainder  of  this  document  is  organized  as  follows.  Chapter  II  provides 
background  on  the  fundamentals  of  incoherent  polarized  light,  the  mathematical 
description  of  its  interaction  with  matter,  and  a  largely  qualitative  discussion  of 
the  practical  aspects  of  polarimetric  remote  sensing.  At  the  end  of  the  chapter, 
there  is  a  brief  discussion  of  atmospheric  effects  on  polarimetric  imaging  systems. 
In  Chapter  III,  the  Cramer-Rao  bound  for  the  estimation  of  polarimetric  imagery  is 
derived  and  applied  in  a  series  of  imaging  case  studies.  The  multichannel  polarimet- 
ric  blind  deconvolution  algorithm  is  derived  and  tested  using  laboratory  data  and 
simulation  in  Chapter  IV.  Next,  the  PCA-based  channel  optimization  procedure  is 
laid  out  in  Chapter  V.  The  document  is  summarized  and  concluded  in  Chapter  VI. 
The  appendices  contain  some  useful  background  derivations  (A),  expansions  on  the 


5 


Cramer- Rao  bound  work  (B),  some  additional  derivation  for  the  blind  deconvolution 
problem  (C),  and  a  mechanical  description  of  the  polarization  imagers  constructed 
for  laboratory  use  (D). 


6 


II.  Overview  of  Polarimetry 


The  purpose  of  this  chapter  is  to  provide  background  on  some  recurring  themes 
throughout  this  research.  Chief  among  these  topics  is  the  mathematical  for¬ 
malism  used  here  to  describe  the  manipulation  of  polarized  light:  the  Stokes  vectors 
and  Mueller  matrices.  These  concepts  are  then  used  to  describe  a  polarimetric  imag¬ 
ing  system.  A  brief  section  is  then  devoted  to  the  phenomenology  of  polarimetric 
remote  sensing.  Finally,  an  overview  of  atmospheric  effects  is  provided. 

2.1  The  Stokes  Parameters 

The  Stokes  parameters  are  a  standard  tool  for  describing  the  polarized  state 
of  light  in  passive  remote  sensing  literature  [6] .  The  utility  of  the  Stokes  parameters 
stems  from  the  fact  that  they  can  be  determined  directly  from  observables.  To  illus¬ 
trate  this  critical  point,  the  Stokes  parameters  are  derived  here  by  closely  following 
the  work  of  Collett  [5]  but  modified  slightly  using  some  ideas  from  statistical  optics. 

The  analytic  signal  representation  of  light  (plane  wave  radiation)  in  a  vacuum 
is  given  by  two  orthogonal  electric  fields,  ux(t)  and  uy(t),  that  are  both  orthogonal 
to  the  direction  of  propagation  [12]: 

ux{t)  =  Ex{t)e~i2™°t  Uy(t)  =  Ey^e-^  (2.1) 

where  Ex(t)  and  Ey(t)  are  complex  valued  functions  that  modulate  the  signal  center 
frequency,  Vo-  This  representation  adequately  describes  light  of  any  bandwidth  and 
any  phase  relationship  including  the  case  where  the  phase  between  the  x  and  y  field 
components  are  partially  or  wholly  dependent  (e.g.  polarized  light). 

The  polarization  state  of  this  field  can  be  determined  by  placing  two  special¬ 
ized  optical  elements  in  its  path  and  measuring  the  result.  The  first  of  these  two 
elements,  a  retarder,  introduces  a  constant  phase  delay,  between  the  x  and  y  field 
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(c)  The  delayed,  attenuated  plane 
wave 

Figure  2.1:  Polarization  transformations  on  a  plane  wave 

components  (figure  2.1b).  The  second  element,  a  polarizer,  transmits  the  portion  of 
the  field  along  angle  6  (measured  with  respect  to  the  x  axis)  and  completely  atten¬ 
uates  the  field  everywhere  else  (figure  2.1c).  The  angle  defined  by  9  is  henceforth 
referred  to  as  the  primary  transmission  axis.  The  resulting  field  along  the  primary 
transmission  axis  is: 

up(t)  —  ux(t)e?^  cos 6  +  uy(t)  sin#  (2.2) 

The  transition  from  (2.1)  to  (2.2)  is  derived  in  Appendix  A.l. 

Since  semiconductor  optical  detectors  respond  to  radiant  intensity,  the  signal 
corresponding  to  equation  (2.2)  at  the  detector  array  is  given  by: 


!{8,<t>)  =  £o  c 


(2.3) 


where  u*  is  the  complex  conjugate  of  up,  (. . .)  represents  the  time  average,  c  is  the 
speed  of  light  and  eo  is  the  permittivity  of  free  space.  (Henceforth,  the  term  eoc  will 
be  suppressed).  Because  v0  is  on  the  order  of  several  terrahertz,  equation  (2.3)  is 
well  approximated  by  an  infinite  time  average.  This  equation  can  be  expanded  into 
its  component  parts: 

J(M)=  (2.4) 

^Ex(t)E*(t)  cos2  6  +  Ey(t)E*(t )  sin2  9 

+  (Ex(t)E*y(t)ej(l>  +  El(t)Ey(t)e-i*)  cosflsin^, 

After  applying  several  well  known  trigonometric  identities: 

/(M)=  (2-5) 

(Ex(t)E*x(t)±±^  +  Ey(t)E*y(t)^^l 

+  (Ex(t)E*(t)  (cos0  +  jsin</>)  +  E*(t)Ey(t)  (cos0  -  j  sin  <j>)) 
and  collecting  like  terms: 

/(M)=  (2.6) 

+  )  (E,(t)E’Jt)  -  Ev(t)E;(t))  cos  (20) 

+  \  (Ex{t)E*(t)  +  ^(^W)  sin  (20)  coscp 
+  |  <£?,(«)£*(«)  -  S-(i)B„(i)>  sin  (20)  sin  4>, 
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the  Stokes  parameters  are  specified  to  be: 


S„  =  {EJt)E-x(t)  +  E,(t)E;(t))  (2.7a) 

51  =  (Ex(t)E-x(t)  -  E,(t)E;(t))  (2.7b) 

52  =  (E,(t)E't(t)  +  E’x(t)E,(t))  (2.7c) 

S3  =  j  {E,(t)E'y(t)  -  El(t)Et(t))  (2.7d) 

such  that  the  signal  at  the  detector  array  is  given  by: 

1(9,  4>)  =  -(So  +  Si  cos  29  +  S2  cos  0  sin  29  +  S3  sin  0  sin  29).  (2.8) 


The  Stokes  parameters  are  therefore  the  intensity  observables  of  the  polariza¬ 
tion  state  of  electro-magnetic  radiation.  Because  it  is  convenient  for  mathematical 
manipulation,  the  Stokes  parameters  are  often  represented  in  vector  notation: 


S 


S0  Si 


S2 


(2.9) 


It  is  important  to  point  out  that  the  Stokes  parameters  were  derived  with 
respect  to  a  specific  x  and  y  coordinate  system.  As  a  result,  agreement  between  the 
raw  measurements  from  a  pair  of  observers  is  wholly  dependent  on  the  orientation 
of  their  reference  frame.  If  the  observers  agree  upon  a  common  coordinate  system 
in  advance  then  their  measurements  can  be  brought  into  agreement  by  applying  the 
proper  transformation  (to  be  described  later  by  equation  (2.17)).  Throughout  this 
work,  the  rc-axis,  or  primary  axis,  will  be  oriented  parallel  to  the  ground  plane  with 
the  y-axis,  or  secondary  axis,  perpendicular  to  it.  In  this  way,  measurements  on  the 
ground  or  in  the  air  will  have  a  common  reference  plane. 
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2.2  Interpretation  of  the  Stokes  Parameters 

By  itself,  the  result  in  equation  (2.8)  provides  little  indication  of  how  the  Stokes 
parameters  should  be  interpreted.  Consequently,  the  purpose  of  this  section  is  to 
introduce  the  relevant  properties  of  the  Stokes  parameters  as  they  pertain  to  this 
research. 

Beyond  the  utility  of  describing  polarization  in  terms  of  intensity,  the  preceding 
derivation  of  the  Stokes  parameters  also  implies  a  number  of  other  important  proper¬ 
ties.  The  use  of  the  analytic  signal  representation  of  polychromatic  electro-magnetic 
radiation  in  equation  (2.1)  demonstrates  that  the  Stokes  parameters  are  equally  well 
defined  for  radiation  of  any  bandwidth.  Consequently,  the  effective  bandwidth  limit 
of  a  Stokes  measurement  is  determined  by  the  experimenter’s  ability  to  impose  the 
same  phase  shift  and  directional  attenuation  at  every  wavelength  in  the  passband. 
The  analytic  signal  representation  also  demonstrates  that  the  polarization  state  of 
the  superposition  of  two  incoherent  fields  is  given  by  the  sum  of  their  Stokes  param¬ 
eters. 

The  physical  interpretation  of  each  of  the  Stokes  parameters  is  also  of  interest. 
Equation  2.7  reveals  Sq  to  be  the  total  intensity  of  the  incident  radiation.  In  other 
words,  Sq  is  the  intensity  that  would  be  observed  if  the  polarizing  elements  were 
removed  from  the  system.  As  a  result,  no  other  Stokes  parameter  can  be  greater 
in  magnitude  than  Sq.  The  remaining  Stokes  parameters  form  a  basis  in  intensity 
space  spanning  all  possible  polarization  states.  Table  2.1  summarizes  the  principle 
direction  for  each  of  the  remaining  Stokes  parameters.  Implied  in  this  table  is  that 
fact  that  Si,  S2,  and  S:i  are  signed  quantities;  S0,  being  the  total  intensity,  is  strictly 
positive. 

In  addition  to  describing  the  orientation  of  polarized  radiation,  the  Stokes 
parameters  also  describe  the  extent  to  which  the  radiation  is  polarized.  Radiation 
with  no  preferred  polarization  state  is  said  to  be  unpolarized.  If  unpolarized  light 
were  to  be  described  using  equation  (2.1),  Ex(t)  and  Ey{t )  would  be  rapidly  varying, 


11 


Table  2.1:  Summary  of  the  Stokes  Parameters 


Normalized 

Stokes  vector  Type  of  polarization 


[i 

[i 

[i 

[i 

[i 

[i 


1  0  o  }T 
-1  0  o  ]T 
0  1  o  }T 
0  -1  o  ]T 
0  0  1  ]T 
0  0  -1 ]T 


linear  horizontal 
linear  vertical 
linear  +45° 
linear  —45° 
right  hand  circular 
left  hand  circular 


independent  random  processes.  This  variation  is  so  fast  that,  during  the  period 
of  any  practical  measurement,  the  instantaneous  polarization  state  of  the  radiation 
would  pass  though  every  possible  polarization  state  many  times.  Hecht  refers  to  this 
state  as  randomly  polarized  or  natural  light  [16]  and  points  out,  along  with  [42], 
that  unpolarized  light  must  be  polychromatic  to  achieve  this  rapid  variation  in  field 
amplitude.  A  consequence  of  this  rapid  amplitude  variation  is  that  the  net  effect  of 
the  retarder  and  polarizer  pair  in  (2.8)  is  the  same  regardless  of  their  orientation. 
By  inspection,  this  condition  requires  that  Si,  S2,  and  S3  be  zero.  Hence  the  Stokes 
vector  of  unpolarized  light  is: 


S 


So 


T 


0  0  0 


(2.10) 


Partially  polarized  light  can  be  considered  to  be  the  sum  of  equation  (2.10) 
with  another  Stokes  vector  that  represents  fully  polarized  radiation.  In  all  cases 
(total,  partial,  or  unpolarized),  the  following  relationship  between  Stokes  parameters 
is  valid: 

S02  >  S2  +  S22  +  S|  (2.11) 
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This  expression  is  bounded  on  two  sides:  total  polarization  exists  when  both  sides 
of  the  expression  are  equal  and  radiation  is  said  to  be  unpolarized  when  the  right 
hand  side  is  zero.  Partial  polarization  occurs  in  every  intermediate  state.  Equation 
(2.11)  leads  to  a  natural  expression  for  the  extent  of  polarization: 

v /s?  +  +  SI 

So 

P  is  referred  to  as  the  degree  of  polarization  and  is  bounded  between  zero  (unpolar¬ 
ized  light)  and  one  (fully  polarized  light).  When  measurements  of  S3  are  unavailable, 
its  value  is  assumed  to  be  zero  and  (2.12)  is  referred  to  as  the  degree  of  linear  po¬ 
larization.  Unlike  the  Stokes  parameters  themselves,  P  is  equivalent  in  all  reference 
frames. 

Additionally,  it  is  also  often  useful  to  express  the  angle  of  linear  polarization, 
ij),  directly  from  the  Stokes  parameters: 

^  =  ^tan“1|^  (2.13) 

which  stems  directly  from  that  fact  that: 


51  =  PS0  cos  2^ 

52  =  PS0  sin  2'ijj 


(2.14) 


Assuming  linearly  polarized  light,  this  last  pair  of  relationships  allows  for  decom¬ 
position  of  /  into  linearly  polarized  and  unpolarized  components  (with  intermediate 
steps  worked  out  in  Appendix  A. 2): 


/  =  -  (1  —  P)  So  +  -  PSq  (1  +  cos  2i/j  cos  29  +  sin  2?/;  sin  29) 
=  +  -V  cos2  (^  -  9) 


(2.15) 
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where  Xu  =  (1  —  P)S0  is  the  unpolarized  component  of  the  signal  and  \p  =  PS0 
is  the  polarized  component;  both  quantities  are  always  greater  than  or  equal  to  0. 
The  idea  of  decomposing  a  source  into  polarized  and  unpolarized  Stokes  vectors  goes 
back  to  Stokes  himself.  A  recent  paper  by  Wolf  explains  the  conditions  under  which 
this  decomposition  is  viable  [49].  Wolf’s  conditions  are  always  met  when  the  target 
is  illuminated  by  natural  light.  In  Chapter  IV,  the  polarized/unpolarized  irradiance 
representation  (rather  than  the  Stokes  vector  representation)  is  employed  to  ensure 
the  statistical  compatibility  of  the  data  model  with  the  measurements. 

2.3  Muller  Matrix  Transformations 

This  section  introduces  the  mathematical  operations  used  to  describe  the  in¬ 
teraction  of  a  Stokes  vector  with  matter.  These  operators  are  called  the  Mueller 
matrices.  Given  a  4  x  4  Mueller  matrix,  T,  the  relationship  between  the  incident 
Stokes  vector,  S,  and  the  transformed  Stokes  vector,  S',  is  given  by  S'  =  T S.  Mueller 
matrices  have  been  defined  for  all  manner  of  optical  elements,  however,  this  discus¬ 
sion  will  be  limited  to  the  the  operators  that  appear  in  this  research:  the  polarizer 
and  the  rotator  [5]. 

A  linear  polarizer,  or  polarization  analyzer,  is  an  optical  element  that  absorbs 
radiation  at  rate  that  is  dependent  on  the  orientation  of  the  incident  electric  field. 
This  held  dependence  is  defined  in  terms  of  a  and  f3  which  are  the  held  transmission 
coefficients  along  the  primary  and  secondary  transmission  axes  as  defined  previously 
in  this  chapter.  Since  a  was  selected  to  be  the  transmission  coefficient  for  radiation 
along  the  primary  axis,  the  following  convention  will  be  adopted  henceforth:  a  >  (3. 
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The  Mueller  matrix  of  a  linear  polarizer,  T,  is  given  by: 

a2  +  j32  a2  —  (32  0  0 

1  a2  —  j32  a2  +  (32  0  0 

T  —  -  (2.16) 

2  0  0  2af3  0 

0  0  0  2a(3 

Polarizers  are  the  central  optical  element  in  this  research  and  are,  therefore, 
worthy  of  a  more  detailed  discussion,  a  and  (3  appear  in  (2.16)  as  second  order  terms 
(e.g  a 2  or  a/3 )  because  the  Stokes  parameters  have  units  of  intensity  which  is  propor¬ 
tional  to  the  square  of  the  mean  electric  field.  Laboratory  quality  polarizers  often 
have  secondary  axis  transmission  coefficients,  /?,  very  near  zero  over  the  bandpass  of 
the  filter.  Conversely,  a  “polarizer”  with  a  —  (3  is  a  neutral  density  filter. 

An  optical  element  that  rotates  the  electric  field  vector  through  an  angle  9  is 
called  a  rotator.  Besides  describing  an  actual  optical  element,  the  rotator  is  also  a 
useful  mathematical  tool  for  translating  Stokes  parameters  between  reference  frames. 
It  is  assumed  that  a  rotator  does  not  attenuate  the  signal  in  any  way  during  rotation 
(i.e.  it  is  an  ideal  rotator).  That  being  said,  the  Mueller  matrix  of  an  attenuating 
rotator  can  be  described  by  the  matrix  product  of  the  ideal  rotator  with  a  polarizer. 
The  ideal  rotation  operator  is  given  by: 

0  0  0 
cos (29)  sin  (26)  0 

(2.17) 

-  sin  (29)  cos (29)  0 
0  0  1 

Thus  the  rotation  operator  also  provides  the  mathematical  means  to  describe  rotation 
of  a  polarizer  by  9: 

T  (9)  =  R-1  (9)TR(9) ,  (2.18) 


1 

0 

R(9)  = 

0 

0 
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which  is  useful  in  practice  because  the  problem  of  characterizing  a  polarizer  in  a  op¬ 
tical  system  can  be  decoupled  into  the  easier  problems  of  measuring  the  transmission 
coefficients  and,  separately,  determining  the  orientation  of  the  device. 

2.4  Mathematical  Treatment  of  Polarization  Imaging 

The  relationship  between  the  at-aperture  Stokes  parameters  and  the  intensity 
sensed  at  each  channel  in  an  imaging  polarimeter  is  described  in  this  section.  This 
derivation  applies  to  any  polarimetric  imager  that  records  each  channel  simultane¬ 
ously.  The  results  provided  here  are  used  extensively  throughout  the  remainder  of 
this  document. 

In  general,  measurement  of  the  Stokes  parameters  requires  four  intensity  mea¬ 
surements  with  9  and  (j)  chosen  such  that  four  unique  instances  of  (2.8)  are  realized. 
In  practice,  the  S3  term  is  approximately  zero  (or,  more  precisely,  (j)  =  0°)  in  remote 
sensing  applications  (see  section  2.5).  Consequently,  only  unique  measurements  are 
required  to  extract  the  relevant  Stokes  parameters  and,  as  shown  below,  the  optical 
element  responsible  for  the  phase  shift  in  (2.2)  may  be  dropped. 

A  channel  model  provides  the  transformation  between  the  incident  Stokes  vec¬ 
tor  image  (or  simply,  Stokes  image),  S,  and  intensity  on  the  detector  array,  repre¬ 
sented  by  Ig.  In  this  arrangement,  6  is  the  orientation  of  the  channel  polarization 
filter  measured  with  respect  to  the  axis  of  the  first  channel.  After  this  transforma¬ 
tion  is  complete,  only  the  total  intensity  reaching  the  detector,  (i.e.  the  Sq  term 
after  transmission  through  the  polarizing  element)  is  of  interest.  Together,  the  total 
intensities  of  each  channel  are  represented  by  an  N  x  1  vector,  I.  Mathematically,  if 

S'(0)  =  T(0)S  (2.19) 
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then 


S'(0)  =  -  [(a2  +  (3 2)  S0  +  (a2  —  /3 2)  cos  (20)  Si  +  (a2  —  (3 2)  sin  (29)  (2.20) 


and 


S'0(  0)  W)  -  S'0(9n_  i) 

lT 

J0  /ei  ...  Ie(N-i) 


(2.21) 


where  the  change  of  variables  from  S'0(9 1)  to  Igi  is  simply  to  avoid  confusion  later 
on  between  the  Stokes  vector  of  the  at-aperture  signal  and  the  actual  intensity  that 
reaches  the  focal  plane.  The  reader  may  compare  this  result  to  the  ideal  case  in 
equation  (2.8). 

Using  equation  (2.20),  a  matrix  transformation,  M,  can  be  defined  to  calculate 
I  directly  from  S: 


M  —  - 
2 


af  +  Pi  (a2  —  (32)  cos  (29i)  (a2  —  /32)  sin  (29i) 


a 


N 


+  Pn  (aN  -  Pn)  cos  (2 On)  (a%  -  (32N)  sin  (2 0N) 


(2.22) 


such  that  I  =  MS.  Note  that  subscripts  have  to  been  added  to  a  and  (3  to  indicate 
that  they  can,  in  general,  vary  from  polarizer  to  polarizer  (though  the  convention 
remains  that  a  >  (3).  Naturally,  the  inverse  (or  pseudoinverse)  of  M  can  also  be 
used  to  calculate  S  directly  from  I: 


S  =  M-1I  (2.23) 

a  and  (3  represent  a  channel  analyzer  in  a  less  than  ideal  case.  That  is  to  say,  the 
preferred  signal  may  be  partially  attenuated  while,  at  the  same  time,  the  remaining 
signal  is  not  completely  suppressed.  As  a  practical  matter,  the  intensity  equation  is 
better  represented  with  weights  on  the  Stokes  parameters  that  are  determined  via 
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lab  calibration  rather  than  separate  measurements  of  a  and  /3: 


M  = 


Q 1  o  •  •  •  a  12 


ClN0  '  '  '  &N2 


(2.24) 


2.5  Polarization  in  Remote  Sensing 

The  purpose  of  this  section  is  to  provide  an  overview  of  polarimetric  remote 
sensing  in  practice.  Besides  providing  a  general  feel  for  the  discipline,  a  primary  goal 
of  this  section  is  to  emphasize  the  point  that  naturally  occurring  objects  with  a  high 
degree  of  polarization  are  unusual.  This  is  the  motivation  for  using  polarization  in 
targeting  and  reconnaissance. 

Polarimetric  imagery  is  highly  dependent  on  the  relative  orientation  of  the 
source,  target,  and  sensor.  Referring  to  figure  2.2,  the  maximum  polarization  re¬ 
sponse  for  natural  materials  tends  to  occur  when  the  phase  angle,  p,  is  near  100°  [7]. 
Likewise,  polarization  is  more  likely  to  occur  in  the  principal  plane  (the  plane  that 
passes  through  the  source  and  target  surface  normal)  than  out  of  it.  Broadening  out 
at  bit,  stronger  polarization  response  is  more  likely  to  occur  in  the  forward  scattering 
direction  than  in  the  back  scattering  (source  side)  direction. 

Plane  (i.e.  linear)  polarization  is  more  likely  to  occur  in  nature  than  elliptical 
polarization  [6,7].  Consequently,  S3  is  approximately  zero  and,  as  such,  is  generally 
not  measured  in  remote  sensing  applications.  To  help  visualize  this  point,  consider 
the  common  case  of  reflection  at  a  smooth  dielectric  surface.  Using  the  familiar 
Fresnel  equations,  Collett  [5]  shows  that  the  Mueller  matrix  transformation  of  this 
physical  process  is  identical  in  form  to  equation  (2.16),  the  Mueller  matrix  of  a  linear 
polarizer.  In  this  case,  a  and  f3  are  determined  by  the  angle  of  incidence  and  the 
indices  of  refraction  of  the  materials  at  their  boundary.  Reflection  on  bare  metals, 
on  the  other  hand,  can  produce  elliptical  polarization  under  the  right  conditions; 
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Surface 


however,  this  and  other  birefringent  materials  are  the  exception  rather  than  the 
rule. 

Multiple  scattering  processes  reduce  the  extent  to  which  light  is  polarized 
[8,16].  Multiple  scattering  processes  in  a  gas  or  aerosol  cloud  can  be  thought  of 
as  the  random  absorption,  re-emission,  and  recombination  of  light.  For  solid  sur¬ 
faces,  especially  those  that  are  semi-transparent,  multiple  scattering  can  be  thought 
of  as  the  net  result  of  many  individual  reflections  from  randomly  oriented  surfaces. 
Consequently,  smooth  surfaces  produce  higher  polarization  states  than  diffuse  sur¬ 
faces  [6].  Scattering  off  of  a  blackbody,  which  is  the  ultimate  diffuse  surface,  has 
no  preferred  polarization  state  even  when  illuminated  with  polarized  light.  Because 
man-made  materials  tend  to  have  smooth,  machined  surfaces,  they  tend  to  reflect 
higher  polarization  states  than  their  surroundings  (water  is  an  obvious  exception  to 
this  rule). 

2.6  Atmospheric  Effects 

The  purpose  of  this  section  is  to  briefly  introduce  the  effects  of  the  atmosphere 
on  polarimetric  imaging  sensors.  Atmospherically  induced  polarization  and  depolar- 
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ization  are  discussed  qualitatively  to  provide  awareness  of  the  issues.  Second,  the 
effects  of  the  atmosphere  on  image  formation  are  provided  as  a  primer  for  the  chapter 
on  blind  deconvolution. 

Rayleigh  scattering  in  the  atmosphere  produces  polarized  sky  shine  [16].  On  a 
molecular  level,  this  scattering  induces  a  preferred  polarization  that  increases  with 
increasing  angle  away  from  the  original  direction  of  propagation.  When  viewed  from 
the  ground,  the  net  effect  of  this  molecular  scattering  is  that  the  sky  takes  on  a 
preferred  polarization  that  is  at  a  maximum  along  azimuths  that  are  perpendicular 
to  the  azimuth  of  the  sun.  Simultaneously,  multiple  scattering  cine  to  atmospheric 
aerosols  and  haze  tend  to  act  as  depolarizing  agents.  In  remote  sensing,  this  polarized 
sky  shine  can  be  reflected  off  of  the  ground  and  into  the  sensor.  The  most  dramatic 
example  of  this  effect  can  be  found  in  the  polarized  reflection  from  ground  targets 
in  shadows. 

In  the  same  way  that  they  effect  sky  shine,  aerosols,  dust,  and  haze  also  tend  to 
diminish  target  polarization  signatures  as  seen  from  an  airborne  sensor.  The  details 
of  this  interaction  are  beyond  the  scope  of  this  dissertation;  however,  the  extent 
of  aerosol  depolarization  is  determined  by  the  type,  depth,  and  number  of  aerosol 
species  in  the  intervening  atmosphere. 

Though  not  a  polarization  effect  per  se,  the  atmosphere  plays  an  important 
role  in  the  formation  of  any  incoherent  image.  In  the  absence  of  atmospheric  effects, 
image  formation  is  given  by  a  convolution  operation: 

i(x)  =  [o  <g>  h]  (x)  (2.25) 

where  h(x)  is  the  sensor  induced  point  spread  function  and  o(x)  is  the  magnified  (or 
minified)  target  image  as  defined  by  geometric  optics.  The  summation  form  of  (2.25) 
can  be  found  at  (4.1).  The  point  spread  function,  even  under  ideal  conditions,  acts  as 
a  low  pass  filter  on  the  geometric  image  and  is  a  strong  function  of  wavelength  and 
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aperture  diameter.  When  the  atmosphere  is  introduced,  uneven  air  temperature 
distributions  cause  spatially  varying  fluctuations  in  index  of  refraction  which,  in 
turn,  modify  the  wavefront  of  the  image  reaching  the  sensor.  The  net  result  of  this 
perturbation  is  a  modified  image  formation  equation: 

i(x)  =  [o®(h  <g)  hatm)\  (. x )  (2.26) 

where  hatm  is  the  atmospheric  operator  on  the  image.  The  grouping  of  the  point 
spread  function  with  the  atmospheric  operator  is  mathematically  unnecessary  but, 
physically,  it  highlights  the  notion  that  the  atmosphere  acts  on  the  image  in  the 
same  manner  that  the  point  spread  function  does. 

Mathematically,  the  atmospherically  corrupted  wavefront  is  represented  as  a 
surface  of  relative  wavefront  phase  delays,  sometimes  referred  to  as  a  phase  screen, 

<p: 


[h  <g>  hatm]  (x) 


fj(p(u)  0— i2i\:kux 


(2.27) 


where  A  is  a  complex  function  that  represents  the  sensor  aperture  and  any  sensor 
dependent  phase  errors.  For  a  diffraction  limited  system,  A{u)  =  1  for  all  u  in  the 
aperture  and  A(u)  =  0  otherwise.  An  explanation  of  what  constitutes  a  typical 
atmospheric  phase  screen  can  be  found  in  [34], 


These  spatial  variations  in  refraction  index  change  over  time.  Consequently, 
hatm  has  a  characteristic  correlation  time,  the  duration  of  which  is  determined  by 
the  turbidity  in  the  atmosphere.  Realizations  of  hatm  are  independent  of  each  other 
from  one  correlation  time  to  the  next.  When  the  exposure  time  of  an  image  is  one 
correlation  time  or  less,  these  index  variations  are  effectively  stationary.  In  this 
regime,  referred  to  as  short  exposure  imaging,  the  dominant  effect  of  the  atmosphere 
is  to  shift  the  image,  though  this  is  also  accompanied  by  other  forms  of  degradation. 
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A  physical  model  of  polarization  imaging  has  been  described  in  this  chapter.  In 
the  following  two  chapters,  this  model  will  be  placed  in  a  statistical  context  through 
which  the  contributions  of  this  research  are  described. 
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III.  Fundamental  Estimation  Bounds  for  Polarimetric  Imagery 


Precise  channel-to-channel  registration  is  a  prerequisite  for  effective  exploitation 
of  passive  polarimetric  imagery.  In  this  chapter,  the  Cramer-Rao  bound  is 
employed  to  determine  the  limits  of  registration  precision  in  the  presence  of  scene 
polarization  diversity,  channel  noise,  and  random  translational  registration  errors 
between  channels.  The  effects  of  misregistration  on  Stokes  image  estimation  are 
also  explored  in  depth.  Algorithm  bias  is  discussed  in  the  context  of  the  bound, 
without  being  estimator  specific.  Finally,  case  studies  are  presented  for  polarization 
insensitive  imagery  (a  special  case)  and  linear  polarization  imaging  systems  with 
three  and  four  channels.  An  optimum  polarization  channel  arrangement  is  proposed 
in  the  context  of  the  bound. 

This  work  was  originally  published  in  [24],  Some  additional  research  on  the 
CRLB  appears  in  Appendix  B.  Included  in  this  appendix  is  a  discussion  of  external 
measurement  versus  joint  estimation,  an  interpretation  of  the  bound  using  corre¬ 
lations,  and  a  method  for  transforming  the  bound  derived  here  so  that  it  may  be 
applied  to  the  model  in  Chapter  IV. 

3.1  Polarimetry  in  the  Context  of  the  Cramer-Rao  Bound 

When  derived  from  imagery,  target  polarization  state  estimates  are  inherently 
prone  to  both  channel  noise  and  registration  errors.  In  this  chapter,  we  examine  the 
combined  effects  of  misregistration  and  channel  noise  in  the  statistical  framework  of 
the  CRLB. 

For  the  present  purpose,  systems  responsive  to  linear  polarization  are  consid¬ 
ered.  Persons  et  al.  noted  that,  among  common  means  of  collecting  polarimetric 
imagery,  translational  and  rotational  registration  errors  between  channels  are  most 
prevalent  [32],  In  this  context,  Persons  and  others  [23,48]  have  made  limited  at¬ 
tempts  to  characterize  the  effects  of  registration  errors  on  polarimetric  imagery  us- 
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ing  established  (but  polarization  insensitive)  registration  algorithms.  In  addition,  a 
recent  registration  study  has  also  been  conducted  for  the  related  problem  of  Mueller 
matrix  imaging  [14].  The  results  of  these  studies  range  from  successful  registration 
within  stated  tolerances  to  complete  failure.  Each  of  these  efforts  is  a  useful  ref¬ 
erence  point,  however,  they  lack  the  depth  required  to  make  complete  comparisons 
across  algorithms  or  imaging  scenarios.  Indeed,  what  is  most  lacking  is  a  theoretical 
context  in  which  to  place  these  results. 

The  first  step  in  providing  this  context  is  to  identify  the  problem  of  estimating 
the  polarimetric  signature  of  a  target  scene  as  a  joint  estimation  of  the  misregistration 
between  channels  and  of  the  scene  itself.  Second,  channel  noise  and  misregistration 
are  random  variables,  therefore,  this  framework  must  be  statistical  in  nature.  Finally, 
this  framework  must  account  for  the  influence  of  deterministic  parameters  such  as 
channel  orientation  and  spacing  as  well  as  polarization  diversity  in  the  scene.  In  this 
paper,  the  CRLB  is  derived  for  this  joint  estimation  problem  and  expounded  as  a 
theoretical  framework  that  meets  all  of  these  criteria. 

The  CRLB  determines  the  lower  limit  on  the  variance  of  any  estimator.  The 
form  of  the  bound  is  different  for  biased  and  unbiased  cases.  With  the  exception 
of  section  3.5,  this  research  focuses  mainly  on  bounds  for  unbiased  estimators.  The 
utility  of  the  CRLB  for  unbiased  estimators  is  well  described  by  Kay  [18];  his  argu¬ 
ments  may  also  be  tailored  to  this  specific  case:  if  a  new  estimator  is  devised,  it  is 
more  efficient  to  compare  it  to  a  single  bound  rather  than  to  attempt  to  simulate 
results  for  many  different  existing  estimators.  Alternatively,  the  bound  can  be  used 
to  determine  if  any  existing  algorithm  meets  this  best  case  scenario.  If  no  such  al¬ 
gorithm  exists,  the  bound  can  motivate  the  search  for  a  better  estimator.  Similar 
studies  via  simulation  (e.g.  Monte  Carlo  methods)  lack  the  theoretical  foundation 
which  the  CRLB  readily  provides.  In  the  present  work,  the  bound  is  also  shown  to 
provide  insight  into  the  relationship  between  the  sensor,  target  scene,  and  estimator. 
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Inspiration  for  applying  this  approach  to  the  polarimetric  imagery  estimation 
problem  comes  from  two  recent  works.  Robinson  and  Milanfar  [33]  applied  the  CRLB 
to  the  problem  of  estimating  the  translational  misregistration  between  pairs  of  (po¬ 
larization  insensitive)  images  in  cases  where  the  underlying  scene  is  not  estimated 
as  a  parameter.  Yetik  and  Nehorai  expanded  this  work  to  describe  limits  for  other 
common  misregistration  issues  (e.g.  rotation,  affine  transformation)  and  for  both 
feature  based  and  intensity  based  registration  algorithms  [50].  Aside  from  insensi¬ 
tivity  to  polarization  effects,  the  underlying  scene  assumption  shared  by  both  sets 
of  authors  obviates  the  need  for  a  joint  estimator  and  is  therefore  inadequate  for  the 
present  purpose. 

The  bound  derived  in  this  work  is  most  applicable  to  so  called  “feature  match¬ 
ing,  area-based”  or  “template  matching”  registration  estimators.  Common  among 
these  methods  is  that  registration  is  achieved  by  comparing  image  intensities  with¬ 
out  regard  to  specific  objects  in  the  scene.  A  survey  of  these  methods  can  be  found 
in  [52], 

The  newly  derived  bound  determines  the  minimum  achievable  variance  on  es¬ 
timates  of  the  translational  misregistration  between  channels  and  of  the  polarization 
state  at  each  point  in  the  scene.  The  inclusion  of  polarization  effects  and  removal 
of  the  known  image  assumption  from  Robinson  and  Milanfar  leads  to  a  bound  cal¬ 
culation  that  is  computationally  impractical  using  direct  methods.  As  such,  matrix 
theory  is  applied  to  the  newly  derived  bound  so  that  extraction  of  the  relevant  terms 
becomes  tenable.  Additionally,  many  well  known  polarization  insensitive  image  reg¬ 
istration  algorithms  exhibit  bias  and,  in  fact,  some  biased  estimators  have  been 
shown  to  outperform  unbiased  estimators  in  theory  and  in  practice.  Though  no  such 
study  has  been  conducted  with  polarimetric  imagery,  a  biased  estimator  form  of  the 
CRLB  is  provided  in  anticipation  of  a  similar  eventual  result. 

The  bound  stands  on  its  own  as  an  algorithm  evaluation  metric,  however,  it 
may  also  be  used  in  a  quantitative  evaluation  of  system  design  criteria.  To  this  end, 
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the  bound  is  used  to  determine  the  conditions  under  which  it  is  necessary  to  generate 
the  channel  registration  parameters  externally  (as  opposed  to  including  them  in  the 
estimator).  Furthermore,  the  polarization  insensitive  case  is  actually  a  special  case 
of  this  new  bound.  As  a  consequence,  Robinson  and  Milanfar’s  bound  is  generalized 
to  the  case  where  the  underlying  image  is  unknown. 

Finally,  the  bound  can  also  be  used  to  suggest  optimal  channel  configurations. 
Several  researchers  have  shown  that  optimum  configuration  (in  terms  of  maximizing 
SNR)  is  achieved  by  distributing  the  channels  evenly  across  all  possible  polariza¬ 
tion  states.  In  the  most  recent  of  these  attempts,  Tyo  improved  on  this  theory  by 
considering  the  optimization  problem  in  the  presence  of  random  channel  calibration 
errors  [46].  In  perhaps  the  most  directly  related  research,  Tyo  also  showed  that 
channels  that  are  evenly  spaced  across  all  linear  polarization  states  are  optimal  in 
a  principal  components  sense  [45].  The  analysis  presented  here  complements  this 
body  of  research  in  that  it  demonstrates  the  same  conclusion  under  different  sets  of 
assumptions  and  using  different  metrics. 

3.1.1  A  Quick  Example.  Building  upon  the  previous  section,  the  results 
of  a  simple  laboratory  experiment  are  presented  in  order  to  illustrate  some  issues 
regarding  the  registration  of  polarimetric  imagery.  The  underlying  assumption  in 
many  registration  algorithms  of  the  aforementioned  “template  matching”  category 
(e.g.  cross-correlation  and  variants)  is  that  the  images  under  test  contain  the  same 
scene  content  but  are  translated  and  corrupted  by  noise.  Other  methods  in  this 
category,  such  as  registration  via  mutual  information,  are  perhaps  more  suited  toward 
the  registration  of  polarimetric  imagery  in  that  a  statistical  dependence  between 
intensity  values  in  each  image  is  assumed  rather  than  a  direct  correspondence  in 
intensity.  To  illustrate  these  differences,  the  following  example  was  contrived  to 
show  how  a  registration  algorithm  can  utterly  fail  when  strong  polarization  content 
violates  its  underlying  assumptions. 
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Figure  3.1  contains  the  result  of  an  attempt  to  register  the  three  channels  of 
an  imaging  polarimeter  using  cross-correlation.  The  target  scene  contains  two  fully 
polarized  bars  from  a  resolution  target  oriented  such  that  the  bars  have  orthogonal 
polarization  states.  The  channels  are  spaced  such  that  each  is  most  sensitive  to 
polarization  at  0°,  60°,  or  —60°.  The  target  polarization  is  rotated  slightly  («  10°) 
out  of  the  sensor  reference  frame  so  that  each  bar  appears  (to  a  greater  or  lesser 
extent)  in  each  channel.  False  color  and  an  arrow  are  used  to  accentuate  the  weak 
signal  of  the  top  bar  in  the  0°  channel. 


(a)  The  misregistered  0°  (b)  The  60°  channel  (c)  The  —60°  channel 

channel 

Figure  3.1:  The  misregistered  polarimetric  images. 

Clearly,  these  three  channels  do  not  contain  the  same  scene  content  which 
is  in  violation  of  the  assumptions  built  into  the  cross-correlation  algorithm.  As  a 
consequence,  the  bottom  bar  in  channel  0°  has  been  misregistered  to  the  top  bars  in 
the  60°  and  —60°  channels.  The  Stokes  parameter  images  of  this  scene  (figure  3.2) 
provide  an  explanation  of  this  behavior. 

Recall  that  Si  and  S'2  can  take  on  negative  values;  the  reader  should  interpret 
the  dark  regions  in  the  Si  and  S'2  images  in  this  way.  This  behavior  is  unlike  So, 
which  is  strictly  positive.  The  comparatively  higher  contrast  in  the  Si  image  occurs 
because  the  signal  in  S2  is  much  weaker  than  in  Si.  If  the  sensor  were  not  sensitive 
to  polarization  then  the  captured  image  in  each  channel  would  be  the  S'o  image. 

Comparing  the  Stokes  images  in  figure  3.2  with  the  intensity  equation  (2.4), 
the  bottom  bar  in  the  0°  channel  image  is  bright  because  its  preferred  polarization 
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(b)  Bars  in  Si 


(a)  Bars  in  So  (b)  Bars  in  Si  (c)  Bars  in  S 2 

Figure  3.2:  The  Stokes  parameter  images  of  the  bar  target. 


state  is  nearly  parallel  to  the  preferred  polarization  state  of  the  channel.  I11  this 
region  I  m  |(So  +  Si)  where  Si  is  a  positive  quantity.  The  small  contribution  from 
S2  in  this  channel  is  ignored.  The  top  bar  is  orthogonally  polarized  to  channel  0° 
and  in  this  region  the  form  of  the  intensity  equation  is  the  same  but  Si  is  negative. 
A  similar  analysis  could  be  conducted  for  the  remaining  channels  but  the  point  of 
this  section  has  already  been  made:  the  phenomenology  of  polarization  imagery  is 
different  than  that  of  traditional  intensity  imagery  and  as  such,  the  rules  developed 
for  image  registration  must  be  reevaluated  in  this  new  context. 


3.2  Bound  Definition  and  Data  Model 

I11  this  section,  the  CRLB  is  defined  and  the  components  required  to  calculate 
the  bound  are  identified.  These  components  are  then  specified  for  the  specific  prob¬ 
lem  of  calculating  the  bound  for  a  joint  registration  and  Stokes  parameter  estimator. 

3.2.1  Definition  of  the  Cramer- Rao  Bound.  Let  Z  be  a  vector  of  random 
variables  that  is  parameterized  by  vector  6.  Define  9  to  be  any  unbiased  estimate 
of  these  parameters  and  z  to  be  one  realization  of  Z.  The  Cramer- Rao  inequality 
provides  the  lower  bound  on  this  estimator’s  error  covariance  matrix  in  terms  of  the 
Fisher  Information  Matrix  (FIM),  J: 


Cov[9 }  >  J-\ 


(3.1) 


where 

j=-£[I>G^(M)t]  (32) 

and  L(9,  z)  is  the  data  log-likelihood  function  [36].  Note  that  E[. . .]  represents  the 
expected  value  operation  over  Z.  The  minimum  variance  for  the  estimate  of  each 
parameter  in  6  is  given  by  the  diagonal  elements  of  J~l.  When  J  is  not  positive 
definite  the  CRLB  is  not  defined. 

The  concept  of  a  log-likelihood  function  may  require  some  elaboration.  A  likeli¬ 
hood  function  describes  how  the  probability  density  for  a  given  measurement  changes 
as  6  changes.  Hence,  calculation  of  the  likelihood  function  requires  knowledge  of  the 
probability  density  function  for  Z,  pg  ( z),  at  each  measured  z.  The  log-likelihood 
function  is  simply  the  natural  log  of  the  likelihood  function 

L(9,z)  =  lnpe  (z) .  (3.3) 

3.2.2  Data  Model.  The  bound  described  in  the  previous  section  is  defined 
in  terms  of  a  set  of  random  variables  Z,  parameter  vector,  9,  and  their  corresponding 
log-likelihood  function,  L(9,  z).  The  purpose  of  this  section  is  to  define  these  items 
for  the  specific  problem  of  generating  polarimetric  imagery  from  noisy,  misregistered 
data. 

In  this  scenario,  the  collected  images  are  realizations  of  Z.  Consistent  with  [33] 
and  [50],  the  images  to  be  registered  are  modeled  as  sampled,  noisy  versions  of 
the  continuously  varying  underlying  target  scene.  Channel-to-channel  point  spread 
function  variations  are  assumed  to  be  minimal,  sampling  is  at  the  Nyquist  frequency 
or  better,  and  noise  in  the  scene  is  zero-mean,  Gaussian  and  IID.  The  sampled 
coordinates  in  the  (arbitrarily  selected)  first  imaging  channel,  f±,  are  taken  to  be  the 
reference  by  which  the  remaining  channels,  f2  to  /at,  are  specified.  Mathematically, 
the  collected  image  z%  for  channel  i  at  pixel  mn  is  given  by: 
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Zi{  m„)  =  /'/ini,,!  + 


(3.4) 


where  mra  is  the  2-dimensional  coordinates  of  a  given  pixel  in  the  image  plane  and 
£j(mn)  is  a  realization  of  the  channel  noise.  The  boundaries  of  mn  are  determined 
by  the  region  over  which  the  collected  images  overlap.  This  overlap  is  assumed  to  be 
square  (1  <  n  <  p2  for  a  p  x  p  image).  In  a  departure  from  the  cited  research,  the 
image  content  in  each  channel  is  determined  by  a  shared  Stokes  parameter  mapping 
rather  than  a  common  intensity  mapping: 


/l(mn) 


2 


E  dijSjinin) 

3= 0 
2 


(3.5) 


fi( mn)  =  E  aijSjtmn  -  v*), 

3=0 

where  v,  is  the  2-dimensional  translational  misregistration  between  fj  and  fi  and  the 
parameter  weights  atj  come  from  equation  2.24. 


For  each  channel,  there  is  a  one-to-one  correspondence  between  the  2-dimensional 
image  and  a  1-dimensional  vector: 


z i  =  [zj(mi) . . .  Zi(mp2)]T 

f i  =  [/.(mi)  •  •  •  fi{ mp2)]r, 

which,  in  turn,  allows  for  a  very  compact  expression  of  the  following  results. 

As  previously  stated,  the  per-pixcl  channel  noise  is  zero  mean,  IID  and  Gaus¬ 
sian;  consequently,  the  data  log-likelihood  function  is  given  by: 

-1  N 

L (0,  z)  =  ^2  (z?:  _  f?;)T  ~  f*)  +  ^  (3-7) 

i=  1 

where  a2  is  the  noise  variance  and  £  is  a  constant  term  that  is  not  dependent  on 
6.  It  is  clear  from  this  equation  that  the  log-likelihood  is  dependent  on  the  region 
of  overlap  defined  by  the  intersection  of  all  images  fj.  This  region  of  intersection  is, 


(3.6) 
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in  turn,  dependent  on  the  relative  misregistration  between  the  images.  Following 
the  lead  of  Yetik,  this  overlap  region  is  assumed  to  be  constant.  The  efficacy  of  this 
assumption  is  greatest  when  the  relative  misregistration  between  images  is  small 
when  compared  to  the  dimensions  of  the  overlap  region;  it  is  reasonable  to  assume 
that  a  multi-channel  polarimeter  will  be  operating  in  this  regime. 


Additionally,  calculation  of  the  Fisher  information  matrix  requires  an  unam¬ 
biguous  ordering  in  the  parameter  vector  9.  Since  the  goal  is  to  place  a  bound  on 
a  joint  estimator  of  the  translational  shifts  between  images  and  the  values  of  the 
Stokes  parameters  at  each  pixel  in  the  image,  this  parameter  vector  will  be  very 
large  indeed: 


9  = 


(3.8) 


where,  analogous  to  the  vectorized  form  of  the  channels  mean  images,  f),  each  Stokes 
parameter  vector  is  given  by: 


Si  =  [Si(mi)...Si(mp2)]T. 


(3.9) 


Scharf  [36]  provides  an  expression  for  the  Fisher  information  matrix  for  m 
observations  of  an  multivariate,  normally  distributed  random  vector  with  covariance 
matrix,  R ,  and  mean,  f: 


_  diT  J  at 

?;R  06: 


(3.10) 


As  shown  in  (3.7),  there  is  one  observation  (i.e.  m  =  1)  for  each  of  N  random  images. 
Each  image  i  has  a  covariance  matrix  R  =  a2 1  and  mean  f It  is  straightforward  to 
show  that  the  FIM  in  this  case  is  simply  a  sum  of  N  instances  of  (3.10): 


J; 


v 


1  <9f(f  dfn 

a2  ^  89 i  89  j 

n= 1  J 


(3.11) 
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3.3  Fisher  Information  for  the  Joint  Estimator 


The  Fisher  information  matrix  from  (3.11)  is  developed  in  this  section  by 
partitioning  J  into  submatrices.  Matrix  partitioning  is  used  to  exploit  the  inherent 
symmetry  in  the  Fisher  information  matrix  and  to  lend  insight  into  the  physical 
interpretation  of  J.  In  [35],  Scharf  and  McWhorter  show  that  a  FIM  in  the  form  of 
(3.11)  may  be  partitioned  so  that  the  subset  FIMs  of  the  parameter  space  may  be 
considered  individually.  In  this  work,  it  is  useful  to  partition  the  parameter  space 
such  the  sub-FIM  V  describes  covariances  amongst  the  translational  registration 
parameters  and  the  sub-FIM  S  describes  covariances  between  Stokes  parameters. 
Define: 


J  = 


V  HT 
H  S 


(3.12) 


where  H ,  to  use  Scharf  and  McWhorter’s  parlance,  relates  the  intercorrelations  be¬ 
tween  the  registration  and  Stokes  partitions  in  9.  In  what  follows,  each  of  V,  H ,  and 
S  are  described  in  detail. 


The  first  matrix,  V,  is  actually  the  FIM  for  an  unbiased  estimator  of  the 
misregistration  between  channels  when  the  underlying  Stokes  images  are  known  a 
priori.  The  inverse  of  V  is  the  CRLB  for  the  registration  estimator  under  this  “known 
prior”  condition.  In  the  two  channel  polarization  insensitive  case,  V  is  the  bound 
from  Robinson  and  Milanfar.  V  is  composed  of  ( N  —  l)2  submatrices  of  the  form: 


Va  = 


J_  8  f  T 

-2  dvi  Li 


0 


2x2 


8  f  T 

a, .  -L I 


if  i  =  j, 
if  i  ±  j 


(3.13) 
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where  02x2  is  a  2  x  2  zero  matrix.  Though  it  does  not  appear  explicitly  in  (3.12).  it 
is  also  useful  to  define  V  with  the  same  form  as  V  but  with  entries: 


Vi 


v 


if  i  =  j, 
if  i  ±  j 


(3.14) 


The  role  of  V  will  be  demonstrated  in  the  following  section.  For  now,  consider  that 
the  form  of  V  demonstrates  that  errors  in  the  estimates  of  the  misregistration  param¬ 
eters  are  uncorrelated  between  channels  when  perfect  knowledge  of  the  underlying 
scene  exists.  When  this  perfect  knowledge  does  not  exist,  V  will  be  used  to  describe 
the  nature  of  the  correlation  between  channels. 


In  the  lower  right  quadrant  of  (3.12)  is  the  3 p2  x  3 p2  matrix  S.  Note  that  S  is 
the  FIM  that  would  be  used  to  estimate  the  bound  on  a  unbiased  Stokes  estimator  if 
the  relative  shifts  between  the  collected  images  were  known  a  priori.  S  divides  into 
3x3  submatrices  of  size  p 2  x  p 2  corresponding  to  combinations  of  Stokes  parameters 
(images) . 


N 


Sjk  9  ( -fp2  X  p-  )  ^  ^  O'ijO'ik 


(3.15) 


i=  1 


where  Ip 2xp2  is  an  identity  matrix  of  rank  p2.  The  block  symmetry  in  S  allows  for 
substantial  further  simplification  by  defining  the  matrix: 


C  = 


Coo  •  •  •  C02 


C20  '  '  '  C22 


N 

j  Cjk  /W  dij^ik 
i— 1 


in  which  case: 


S  ^2  C  Ip2xp2 


(3.16) 


(3.17) 
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where  (8)  is  the  Kronecker  product  [13].  This  Kronecker  product  representation 
significantly  simplifies  inversion  of  the  S  matrix. 


S  1  =  (J2C  1  (8)  ip2xp2 


(3.18) 


In  other  words,  the  inversion  of  the  3 p2  x  3 p2  matrix  S  is  solved  by  simply  inverting 
the  3x3  matrix  C. 


Connecting  the  matrices  V  and  S  is  the  3 p2  x  2 (IV  —  1)  matrix  H .  Physically, 
if  H  were  the  zero  matrix  then  the  bounds  on  the  registration  and  Stokes  param¬ 
eters  could  be  determined  independently  of  each  other.  Proof  of  this  statement  is 
provided  in  section  (3.4).  In  the  process  of  defining  H,  it  becomes  obvious  that 
this  independence  condition  is  never  met.  H  is  composed  of  3  x  (N  —  1)  readily 
identifiable  submatrices: 


TT  _  aji 
Hlj  ~  U2 


d 


f/ 


dvj  J 


(3.19) 


Unless  the  collected  image  is  constant  everywhere  (i.e.  the  derivative  of  the  image  is 
zero  everywhere)  then  H  is  non- zero  and  the  covariance  bounds  must  be  determined 
jointly.  More  on  the  relationship  between  H,  V  and  V  can  be  found  in  Appendix 


B.l. 


3-4  Bound  Derivation 

The  Fisher  information  matrix  has  been  shown  to  contain  [2(1V  —  1)  +  3 p2]  x 
[2(1V  —  1)  +  3p2]  entries.  To  put  the  enormity  of  this  matrix  into  perspective,  consider 
the  bound  calculation  for  a  four  channel  polarimeter  used  to  estimate  the  first  three 
Stokes  parameters.  Assuming  a  512x512  overlap  region,  the  corresponding  FIM  has 
approximately  6.18  x  1011  entries.  Inversion  of  such  a  large  matrix  is  prohibitive. 
In  this  section,  the  partitioning  of  the  Fisher  information  matrix  from  the  previous 
section  is  exploited  to  make  this  inversion  problem  tractable. 
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3-4-1  Block  Matrix  Inversion.  As  previously  discussed,  the  variance  bound 
for  each  parameter  in  9  is  given  by  the  diagonal  entries  of  J~1.  Consequently, 
computational  expense  can  be  significantly  reduced  by  avoiding  the  unnecessary 
calculation  of  many  of  the  off  diagonal  terms  in  the  inverse.  A  trivial  rearrangement 
of  the  partitioned  inverse  of  a  block  matrix  in  [15]  provides  the  following: 

Bv=(V-  HTS~lH ) (3.20) 

and 

Bs  =  S~l  +  A-1  ( HBvHt )  S -1  (3.21) 

where  Bv  and  Bs,  which  are  submatrices  of  J~1,  identify  the  bounding  covariance 
matrices  on  the  shift  parameters  and  Stokes  images.  As  an  aside,  if  H  =  0  then 
Bv  =  l/-1  and  Bs  =  S'-1,  thus  demonstrating  the  physical  interpretation  of  H  put 
forth  in  the  previous  section. 

Equations  (3.20)  and  (3.21)  readily  demonstrate  how  uncertainty  in  the  mis¬ 
registration  between  channels  impacts  uncertainty  in  the  Stokes  image  estimates  and 
vice  versa.  Similar  to  the  H  —  0  case,  if  the  underlying  Stokes  images  are  known  per¬ 
fectly  then  the  CRLB  on  the  misregistration  estimates  would  simply  be  Bv  =  V -1. 
Likewise,  if  perfect  knowledge  of  the  registration  parameters  existed  then  B$  =  S'-1. 
Consequently,  it  is  clear  that  Bv  and  Bs  are  always  larger  than  V~x  and  S in  the 
absence  of  perfect  knowledge. 

3-4-2  Simplified  Registration  Parameter  Bound.  Matrix  Bv  is  addressed 
first  because  it  is  required  to  calculate  Bs-  As  preliminary  work,  note  that: 

HTS~lH  =  a2 Ht  {C~l  <g>  Ip 2xp2)  H.  (3.22) 
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Consequently,  HTS  1H  can  itself  be  partitioned  into  a  matrix,  D,  such  that: 

2  2 

<3’23) 

i=0  j= 0 

which  provides  the  opportunity  to  profitably  apply  equations  (B.l)  and  (B.2)  from 
Appendix  B.l: 

E  E  Clj1  ahiahjVhh  if  h  =  k, 

Dhk  =  i°j=°  _  (3.24) 

E  E  C^ahiakjVhk  if  h^k 

2=0  j= 0 

Furthermore,  it  is  mathematically  expedient  to  add  and  then  subtract  V  into  the 
definition  of  D  such  that: 

D  =  —Wv  •  (v  +  y)  +  V  (3.25) 

where  •  is  the  Hadamard  product  and  all  sensor  dependent  terms  have  been  bundled 
into: 

BA  =  (l(jv— i)x(iv— i)  —  BC  TAT )  <8)  12X2  (3.26) 

The  matrix  A  is  formed  of  the  2  to  N  rows  of  M ,  the  matrix  of  per  channel  Stokes 
weighting  parameters  from  (2.24).  The  bound  on  the  shift  estimates  can  now  be 
expressed  concisely: 

Bv  =  {V-  D)~l  =  Wv  •  (v  +  y)  1  (3.27) 

As  anticipated,  V  has  indeed  provided  the  constituents  for  the  cross-covariance  terms 
in  Bv. 

3-4-3  Simplified  Stokes  Parameter  Bound.  The  impetus  behind  much  of  the 
preceding  section  was  to  avoid  having  to  manipulate  an  unwieldy  Fisher  information 
matrix  directly  while  still  achieving  the  CRLB  on  the  registration  parameters.  The 
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goal  in  this  section  is  to  do  the  same  for  the  bound  on  the  Stokes  image  estimators, 
Bs ■  We  propose  that  the  variance  bound  on  the  Stokes  parameter  estimate  for  any 
one  pixel  in  the  image  is  of  less  interest  than  the  average  bound  across  the  image. 
In  turn,  this  calculation  is  significantly  simplified  by  applying  properties  of  the  trace 
and  of  the  Kronecker  product  as  shown  in  Appendix  B.2,  where  the  bound  is  derived 
in  detail.  What  follows  are  the  highlights  of  this  derivation. 

Define  B$i,  the  covariance  matrix  for  an  estimator  of  Stokes  image  S;,  to  be 
a  submatrix  of  Bg.  Equivalently,  let  Tj  be  the  submatrix  of  S  corresponding  to  the 
Stokes  image  Sy 

IT1  =  ^2C^®W  (3-28) 

and  define: 

=  cr2^"1  ®  Ip2xp2  (3.29) 


where 


n- 1 

°i0 


then  the  average  bound  for  Stokes  image  S *  is  defined  to  be: 


(3.30) 


(BSi)  =  \tr  (. BSi ) 
jr 

Expanded  out,  the  trace  term  is: 


(3.31) 


tr  (. BSi )  =  tr  (rr1)  +  tr  [d*-1  (. HBVHT )  $“T; 


(3.32) 


and  because  of  the  implicit  symmetry  of  V  and  V : 


tr  [&-1  ( HBvHt )  $“T]  =  a2vec 


=  a  tr 


Ws,;  •  (V  +  V 

WSt  •  (v  +  v)  B 


vec  (Bv 


(3.33) 
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where 


BA*  =  A  (C^Cr1)  AT  ®  l2x2  (3.34) 

this  result  can  be  folded  back  into  (3.31)  to  produce: 

2 

(BSi)  =  a2C^  +  °^tr  [wSi  •  (v  +  y)  £?,]  (3.35) 

Hence,  substantial  simplification  of  the  bound  calculation  is  achieved.  The  similar¬ 
ities  between  (Bsi)  and  Bv  are  clear.  The  terms  in  B$  that  are  not  contained  in 
Bgi  can  be  ignored  because  they  do  not  influence  the  trace.  Also,  it  is  interesting 
to  note  how  the  sensor  itself  (realized  by  HA*  and  Wv  in  Bv)  plays  opposing  roles  in 
equation  (3.35)  analogous  to  multiplication  and  division  if  this  were  a  purely  scalar 
case. 

3.5  Biased  Estimators 

A  study  of  estimator  bias  in  the  registration  of  polarimctric  imagery  has  not 
been  addressed  in  the  literature.  However,  there  are  several  well  known  examples 
of  estimator  bias  in  traditional,  polarization  insensitive  registration  algorithms.  In 
anticipation  of  this  future  work,  what  follows  is  a  discussion  of  incorporating  bias 
into  the  results  from  sections  3.4.2  and  3.4.3. 

According  to  Scharf  [36],  the  CRLB  on  the  covariance  matrix  of  a  biased  esti¬ 
mator  is  given  by: 

Cov  [fll  >  AtJ"1A  (3.36) 

where,  J,  is  the  Fisher  information  matrix  (identical  to  the  unbiased  case),  6  is  the 
estimate  of  6  and 

A=dB[C  (3’37) 

is  the  partial  derivative  of  the  expected  value  of  the  estimator.  The  expected  value 
of  6  is  estimator  specific  (i.e.  registration  algorithm  specific)  and,  furthermore, 
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the  CRLB  of  a  biased  estimator  may  be  higher  or  lower  than  that  of  an  unbiased 
estimator.  Note  that  in  the  unbiased  case,  A  =  I  and  equation  (3.36)  reduces  to 
equation  (3.1).  In  what  follows,  a  specific  partition  of  A  is  defined  via  subscript. 
For  instance,  A„  refers  to  a  partition  of  A  dealing  strictly  with  estimates  of  the 
registration  parameters  whereas  As0  refers  to  estimates  of  the  Stokes  image  Sq. 

The  biased  estimator  bound  for  the  registration  parameters,  Bv ,  is  easily 
achieved  by  combining  (3.27)  with  (3.36): 


b„  =  a: 


wv 


V  +  v 


-1 


A, 


(3.38) 


which  is  made  possible  by  exchanging  J _1  in  (3.36)  with  its  submatrix  Bv  (from  the 
unbiased  case).  In  what  appears  to  be  a  very  small  step,  this  equation  shows  that 
the  CRLB  can  be  decomposed  into  an  scene  specific  part,  V  +  V,  a  sensor  specific 
part  Wv,  and  an  estimator  specific  part,  A„.  This  separation  may  not  be  complete 
in  that  the  channel  spacing  effects  V  +  V  and,  more  than  likely,  bias  in  the  estimator 
will  be  to  some  extent  affected  by  scene  polar imetric  content. 

Without  the  context  of  a  specific  registration  algorithm,  interpretation  of  the 
biased  CRLB  for  a  Stokes  estimator  is  less  straightforward.  Combining  (3.32)  with 
(3.36)  reveals  that  the  trace  on  the  biased  Stokes  bound  is: 


tr  (6a, )  =  tr  (A&TpA*  +  Aj.t"1  (HBvHr)  (3.39) 

which  may  be  reduced  somewhat  by  substituting  in  for  <h“1  in  Appendix 

B.2  and  noting  that: 


<!>„•  T AsiAl^fe;  1  =  a4C,  TC 1  (S>  A.<?,;Ac, 


(3.40) 
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resulting  in: 


tr  (£Si)  =  tr{ ATSiF^ASi)  +  aHr  ( HT  {C~TC^  ®  #£„)  (3.41) 

which  is,  if  not  quite  as  simple  as  the  registration  parameter  case,  certainly  more 
computationally  efficient  than  directly  inverting  the  FIM.  In  terms  of  interpretation, 
the  trace  is  composed  of  the  sum  of  two  terms,  the  first  of  which  depends  only  on  the 
estimator  and  the  sensor.  Later,  in  section  3.6.2,  it  is  shown  that  the  equivalent  to 
this  term  in  the  unbiased  case  dominates  the  solution  to  the  problem  of  an  optimum 
channel  spacing.  It  will  be  left  to  future  work  to  show  that  this  is  the  case  for  specific 
biased  estimators. 

Though  there  are  currently  no  studies  of  bias  in  the  registration  of  polarimetric 
imagery,  the  interested  reader  will  be  well  served  by  referring  to  [33]  for  an  in  depth 
discussion  of  bias  in  gradient  based  registration  algorithms  for  polarization  insensitive 
imagery. 

3.6  Example  Bound  Calculations 

In  the  examples  that  follow,  an  unbiased  estimator  is  assumed  throughout. 

3. 6. 1  Bounds  on  Polarization  Insensitive  Imagery.  Bounds  on  polarization 
insensitive  imagery  are  presented  here  as  a  special  case  of  the  results  derived  in 
section  3.4.  This  case  serves  as  both  a  practical  example  of  how  the  CRLB  can  be 
applied  effectively  and  as  a  bridge  between  this  work  and  the  work  of  Robinson  and 
Milanfar. 

We  seek  the  bound  on  a  joint  estimator  of  a  scene  and  the  registration  pa¬ 
rameters  amongst  N  noisy  samples  of  this  image.  Recall  that  S(i  represents  total 
scene  intensity,  therefore,  So  is  the  only  Stokes  image  of  interest.  In  this  scenario, 
the  matrix  in  equation  (3.16)  reduces  to  a  scalar,  C  =  N,  because  there  is  only  one 
parameter  per  pixel  to  be  estimated.  The  channel  transmission  coefficient,  al0 ,  is 
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assumed  to  be  unity  for  each  channel.  In  that  case,  M  =  l(N-i)xi-  Therefore,  each 
of  Wv  and  Ws0  are  the  2  (N  —  1)  x  2  (N  —  1)  matrices: 


and 


j1  -  jr  if  *  =  h 

\~x  if'  *'  /  j 


(3.42) 


w*  =  h 


l-2(iV— l)x2(iV— !)• 


(3.43) 


To  reiterate  an  important  point,  both  V  and  V  are  image  dependent.  That 
being  said,  the  overall  bound’s  behavior  with  increasing  N  can  be  predicted  by  the 
behavior  of  Wv  and  Ws0  for  any  image.  Equation  (3.42)  shows  that  Bv  =  2V~1  in 
the  two  channel  case  and  Bv  =  V _1  in  the  limit  of  N.  Essentially,  the  V  matrix 
is  suppressed  by  the  1/N  terms  in  Wv  as  N  increases.  These  suppressed  terms 
represent  the  decreasing  influence  of  each  individual  image  as  the  true  underlying 
image  emerges.  In  this  limit,  the  bound  parameters  for  each  image  asymptotically 
achieve  the  bound  predicted  by  Robinson  and  Milanfar  in  their  2004  work.  Assuming 
a  few  images  do  not  differ  substantially  from  the  rest  of  the  ensemble  (e.g.  due  to  a 
large  translational  error  or  parallax)  then  the  results  follow  a  (l  —  ^)  1  progression. 

Likewise,  the  progression  of  the  estimated  image  toward  the  true  image  can  be 
ascertained  from  the  behavior  of  ( Bs0 )■  Again,  we  consider  the  endpoints.  In  the 
two  channel  case: 

(Bs0)\n=2  =  ^  +  y  (3.44) 

P  Z 

and  the  large  N  case: 


(Bs 0)| 


N^large 


a2  2  (N  —  1)  u2 
jR  W2  +  TV 


(3.45) 
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2 

It  would  appear  that  the  behavior  of  the  this  bound  is  dominated  by  ^ .  As  shown 
in  figure  3.3,  the  broad  applicability  of  the  preceding  observations  can  be  verified 
experimentally. 

Twenty  realizations  of  two  test  images,  each  with  unique  spatial  characteristics, 
are  generated  with  white  Gaussian  noise  statistics  and  an  uniformly  distributed 
random  shift  error  between  ±5  pixels  in  any  direction.  The  scene  average  signal- 
to-noise  ratio  is  3.  The  bound  on  the  average  pixel  value  estimates,  ( Bs0 ),  and 
the  average  misregistration  bound  (i.e.  tr(Bv)/(N  —  1))  are  normalized  and  plotted 
against  the  number  of  frames,  N.  Normalization  is  carried  out  to  illustrate  the  A 
behavior  of  the  intensity  variance  bound  and  the  (l  —  A)  1  behavior  of  the  shift 
estimator  bound.  For  comparison,  each  of  these  predicted  curves  are  shown  as  a  red 
dashed  line  underneath  the  data. 

The  expected  behavior  of  the  bound  with  increasing  N  is  confirmed  for  these 
disparate  examples.  The  deviation  of  average  registration  parameter  bound  in  the 
G.G.  Stokes  image  from  the  expected  trend  demonstrates  the  slight  influence  of 
changes  in  image  overlap  area  on  the  results. 

3.6.2  The  Four  Channel,  Three  Stokes  Case.  Unambiguous  determination 
of  a  linear  polarization  state  requires  three  or  more  polarization  imaging  channels. 
In  the  following  section,  we  demonstrate  that  a  minimum  of  four  channels  is  required 
for  the  joint  registration  and  Stokes  parameter  estimation  problem.  In  this  section, 
example  Cramer-Rao  bounds  are  calculated  as  a  function  of  polarization  channel 
orientation  for  three  distinct  imaging  scenarios.  The  purpose  of  this  exercise  is  to 
illustrate  the  effects  of  channel  orientation  on  the  CRLB  and,  in  this  context,  to 
describe  the  optimum  channel  configuration. 

Along  with  scene  content,  the  orientation  of  each  polarization  channel  can  be 
expected  to  affect  the  bound  of  the  estimated  Stokes  parameters.  As  discussed  in  the 
introduction,  Tyo  has  shown  that  an  optimum  combination  of  polarization  channels 
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Figure  3.3:  Polarization  insensitive  imagery  and  its  estimation  bounds 


exists  in  a  principal  components  sense  for  monochromatic,  fully  polarized  radiation 
that  is  uniformly  distributed  over  all  states  of  linear  polarization  (e.g.  not  image 
specific)  [45].  Tyo’s  optimum  combination  of  channels,  in  the  four  channel  case, 
corresponds  to  a  spacing  of  45°.  Since  his  work  was  brought  about  from  different 
assumptions  and  desired  outcomes,  it  is  of  interest  to  compare  this  result  to  the 
optimum  configuration  as  defined  by  the  CRLB. 

Examples  from  three  real-world  polarization  imagers  are  selected  to  represent 
the  diversity  of  polarization  imaging  scenarios.  The  first  case  is  the  bar  target  from 
figure  3.2  which  was  collected  in  the  lab  using  a  three  channel  imaging  pola.rimeter 
in  the  visible  regime.  The  ultraviolet  astronomical  data  of  galaxy  Markarian  3  comes 
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from  the  (now  decommissioned)  Faint  Object  Camera  that  was  aboard  the  Hubble 
Space  telescope  until  2002.  Finally,  the  machine  vision  example  of  a  printed  circuit 
board  is  taken  from  laboratory  data  collected  in  the  visible  regime  (using  a  different 
polarimeter  than  in  the  the  bar  target  case).  The  Stokes  images  for  the  Markarian 
3  and  PC  board  cases  are  shown  in  the  following  figures. 


(a)  Markarian  3  in  Sq  (b)  Markarian  3  in  Si 

Figure  3.4:  The  Markarian  3  test 


scene 


(a)  PC  board  in  So  (b)  PC  board  in  Si  (c)  PC  board  in  S2 
Figure  3.5:  The  printed  circuit  board  test  scene 


The  bound  characteristics  are  determined  by  varying  the  angular  spacing  be¬ 
tween  polarization  elements  in  one  degree  increments.  The  separation  angle  between 
channels  is  measured  with  respect  to  the  first  channel,  which  is  held  fixed.  Each 
channel  is  evenly  spaced  and,  in  the  plots  that  follow,  this  channel  spacing  is  used 
to  index  the  bound  results.  For  example,  the  reader  can  infer  that  a  channel  spacing 
of  5°  corresponds  to  an  orientation  of  5°,  10°,  and  15°  for  the  second,  third,  and 
fourth  channels  with  respect  to  the  first.  Figure  (3.6)  shows  the  average  bound  on 
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the  Stokes  parameters  for  each  of  the  three  test  cases.  Each  plot  is  normalized  by 
o~2  since  it  may  be  divided  out  of  equation  (3.35). 

As  should  be  expected,  the  largest  variances  in  Stokes  image  estimates  occur 
when  the  polarization  channels  are  closely  spaced.  Somewhat  less  expected  is  the 
observation  that  scene  content  appears  to  have  little  influence  on  individual  results 
for  the  Stokes  estimator  in  terms  of  the  overall  trends.  Though  there  are  small  scale 
differences,  the  Stokes  parameter  bounds  for  all  images  closely  follow  the  evolution  of 
CT1  with  increasing  channel  spacing.  This  result  is  significant  for  two  reasons.  First, 
there  appears  to  be  global  agreement  as  to  which  channel  spacing  is  most  preferable, 
at  least  for  the  test  cases  sampled  here.  Second,  it  would  appear  that  each  Bs%  is 
very  well  approximated  by  direct  interrogation  of  S'-1,  at  least  in  an  average  sense. 
In  other  words,  the  Stokes  parameter  bounds  are  effectively  scene  independent  with 
respect  to  channel  spacing.  Recall  that  S'-1  is,  by  itself,  the  CRLB  on  the  Stokes 
estimates  when  perfect  knowledge  of  the  misregistration  parameters  is  available, 
consequently,  these  observations  apply  equally  to  that  scenario. 

The  bound  for  the  S \  and  S2  estimates  meet  at  Tyo’s  predicted  optimal  channel 
spacing  of  45°,  however,  Si  has  a  minimum  bound  at  a  somewhat  closer  channel 
spacing.  Consequently,  there  is  no  global  minimum  bound  for  all  Stokes  parameters. 
Rather,  the  45°  spacing  is  the  point  where  there  is  no  preferred  parameter.  This 
point  is  significant  because,  as  stated  previously,  the  Stokes  parameters  are  defined 
with  respect  to  some  coordinate  system  and,  as  this  system  changes  in  relation  to 
the  target  (e.g.  through  camera  motion)  then  scene  content  can  shift  between  S± 
and  S2.  With  this  qualifier,  a  45°  channel  spacing  can  be  said  to  be  optimal  for  a 
four  channel  system. 

Figure  (3.6)  also  contains  a  plot  of  the  average  registration  parameter  bound 
for  each  of  the  three  channels.  As  before,  a 2  has  been  normalized  out.  Unlike  the 
bounds  on  the  Stokes  parameters,  these  bounds  depend  both  on  scene  content  and 
channel  orientation.  Consistent  with  the  Robinson  and  Milanfar  analysis  of  V  in 
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the  polarization  insensitive  two  channel  case,  the  difference  in  bound  magnitude 
correlates  with  the  amount  of  high  spatial  frequency  content  in  the  test  images. 
The  PC  board  image,  with  its  multi-faceted  geometric  features,  generates  the  lowest 
bound  while  the  galaxy  Markarian  3,  which  has  largely  diffuse  features,  generates 
the  highest  bound. 


(a)  Markarian  3  Stokes  bounds 


(c)  PC  board  Stokes  bounds 
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Figure  3.6:  Bound  results  for  the  four  channel  case 


3.6.3  The  Three  Channel,  Three  Stokes  Case.  The  three  channel  case  is 
important  because  it  is  the  minimum  number  of  channels  required  to  completely 
describe  linear  polarization  (i.e.  So,  Si,  S2).  The  interested  reader  may  easily  verify 
this  three  channel  prerequisite  for  themselves  by  attempting  to  calculate  Bs  with 
only  two  channels.  Three  channel  polarimeters  are  also  common  in  practice;  all 
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examples  in  the  previous  section  were  originally  collected  by  different  three  channel 
systems.  In  this  section,  the  three  channel  polarimeter  is  examined  from  a  joint 
estimation  perspective. 

First,  consider  the  bound  matrix  Bv  for  some  combination  of  three  polarimeter 
channels  represented  by  weight  matrix  A.  Immediately,  a  mathematical  difficulty 
arises: 

Wv  =  04x4  (3.46) 

in  which  case  equation  (3.27)  is  uninvertible  or,  in  other  words,  the  CRLB  for  the 
joint  estimator  (biased  or  unbiased)  is  infinity. 

At  first  glance,  the  result  in  (3.46)  appears  to  be  inconsistent  with  the  fact  that, 
three  channel  polarimeters  are  routinely  employed  in  practice.  The  key  difference, 
however,  is  that  the  data  from  these  systems  are  not  reduced  using  joint  estimation 
algorithms.  Rather,  shift  estimation  and  Stokes  estimation  are  treated  as  separate 
problems.  In  other  words,  external  measurement  of  the  registration  parameters  is 
always  required  in  the  three  channel  case. 

In  terrestrial  remote  sensing  applications,  it  is  often  the  case  that  much  of  the 
scene  is  dominated  by  weakly  polarized  content  (an  obvious  exception  is  a  target 
scene  composed  over  the  water).  In  this  case,  the  channel-to-channel  registration 
problem  for  each  pair  of  images  is  well  approximated  by  the  two  channel  polarization 
insensitive  case  described  in  section  (3.6.1). 

Assuming  external  measurement  of  the  registration  parameters  by  whatever 
means,  there  is  likely  a  regime  where  S'-1  will  dominate  equation  (3.21),  analogous 
to  the  four  channel  joint  estimation  examples.  Following  the  prescribed  method  in 
the  previous  section,  the  diagonal  terms  in  Jjj1  will  be  jointly  minimized  for  S i  and 
S2  when  the  angular  separation  between  channels  is  60°.  In  this  sense,  the  CRLB 
for  the  Stokes  estimates  can  be  compared  across  the  three  and  four  channel  cases. 
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Table  3.1: 
regime 


Average  bounds  on  Stokes  parameter  estimates  in  the  S  1  dominant 
Parameter  3-channel  case  4-channel  case 


50 

51 

52 


L2 

2 


k2 

kl 


Note  from  table  3.6.3  that  So  follows  the  A  pattern  established  for  the  polarization 
insensitive  case  while  Si  and  S2  go  by  A. 


3. 7  Chapter  Summary 

In  this  chapter,  the  Cramer-Rao  lower  bound  is  derived  for  the  problem  of 
jointly  estimating  Stokes  images  and  random  misregistration  parameters  in  the  pres¬ 
ence  of  channel  noise.  Direct  inversion  of  the  prohibitively  large  Fisher  information 
matrix  is  bypassed  by  applying  matrix  theory  to  express  the  resulting  bounds  in  a 
tractable  form.  The  effects  of  estimator  bias  on  the  bound  are  discussed  up  to  the 
point  where  it  becomes  necessary  to  identify  a  specific  estimator.  The  bound  is  then 
used  to  describe  three  and  four  channel  polarimetric  imaging  systems  as  well  as  the 
special  case  of  registering  N  polarization  insensitive  images.  The  bound  itself  is  a 
useful  evaluation  metric  because  it  incorporates  both  sensor  and  estimation  algo¬ 
rithm  effects  and  can  be  used  to  describe  these  effects  theoretically  in  a  way  that 
pure  simulation  can  not.  In  addition,  the  following  general  conclusions  can  be  drawn 
from  the  results: 

•  The  minimum  achievable  estimator  variance  is  scene  and  channel  orientation 
dependent  in  the  unbiased  case.  The  estimator  itself  also  effects  the  bound  in 
the  biased  case. 

•  The  bounds  on  the  registration  and  Stokes  parameters  are  dependent. 
Furthermore,  specific  conclusions  can  be  drawn  from  the  section  3.6  case  studies: 
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•  In  the  polarization  insensitive  case,  the  CRLB  derived  by  Robinson  and  Mi- 
lanfar  [33]  is  shown  to  be  the  asymptotic  limit  of  the  joint  estimation  bound 
as  the  number  of  frames  approach  infinity. 

•  In  the  three  channel  case,  the  CRLB  for  any  joint  estimator  is  infinity.  Con¬ 
sequently,  the  registration  parameters  and  Stokes  images  must  be  estimated 
separately. 

•  The  form  of  the  bound  suggests  that  the  optimum  channel  spacing  is  60°  and 
45°  respectively  in  the  three  and  four  channel  cases.  Optimum,  in  this  context, 
refers  to  a  joint  minimum  bound  for  ,Sj  and  S^.  These  arrangements  do  not 
guarantee  that  the  bound  on  the  registration  estimator  is  minimized. 
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IV.  Blind  Deconvolution  of  Polarimetric  Imagery 


A  maximum  likelihood  blind  deconvolution  algorithm  is  derived  for  incoherent 
polarimetric  imagery  using  expectation  maximization.  In  this  approach,  the 
unpolarized  and  fully  polarized  components  of  the  scene  are  estimated  along  with  the 
corresponding  angles  of  polarization  and  channel  point  spread  functions.  The  scene 
state  of  linear  polarization  is  determined  unambiguously  using  this  parameterization. 
Results  are  demonstrated  using  laboratory  data  and  simulation. 

The  preponderance  of  this  chapter  was  originally  published  in  [25].  In  an 
expansion  of  this  original  work,  sections  4.6  and  4.7  include  an  algorithm  comparison 
and  an  attempt  to  better  understand  polarization  angle  estimation,  both  in  the 
context  of  simulation.  To  bridge  this  chapter  with  the  last,  the  Cramer- Rao  bound 
for  the  ML  estimator  derived  here  can  be  found  in  the  appendix,  section  B.5. 

4-1  Polarimetric  Image  Restoration 

However  polarization  information  is  conveyed,  it  is  generally  the  result  of  a 
fusion  of  multiple  latent  images,  each  collected  through  a  different  polarization  an¬ 
alyzer.  In  the  presence  of  strong  polarization  features,  these  channels  may  generate 
images  with  substantial  content  variations  between  them.  In  this  sense,  the  fusion 
of  polarization  imagery  is  similar  to  fusion  across  spectral  regimes  or  perhaps  among 
medical  imaging  modalities.  To  some  extent,  these  analogies  are  imperfect  because 
the  relationship  between  channels  in  polarimetric  imagery  can  be  accurately  modeled 
without  knowledge  of  the  target  itself. 

Like  all  imaging  systems,  the  sensor  limits  the  ability  of  the  image  to  faithfully 
reproduce  the  intended  target.  These  limits  are  manifested  by  noise  and  by  the  blur¬ 
ring  effects  of  the  channel  point  spread  function.  Additionally,  it  is  often  the  case 
in  remote  sensing  that  the  response  of  each  channel  varies  due  to  the  wavefront  cor¬ 
rupting  effects  of  the  atmosphere.  The  problem  of  estimating  the  true  target  image 
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in  the  presence  of  unknown  atmospheric  corruption  is  known  as  blind  deconvolution 
and  has  been  studied  extensively  in  the  traditional  imaging  (i.e.  polarization  insensi¬ 
tive)  case.  An  excellent  overview  of  the  discipline  can  be  found  in  references  [21,22]; 
however,  there  are  also  more  recent  examples  in  the  literature.  In  this  paper  we 
derive  and  test  a  maximum  likelihood  blind  deconvolution  algorithm  for  incoherent 
polarimetric  imagery. 

The  proposed  blind  deconvolution  algorithm  is  a  generalization  of  the  work  of 
Schulz  on  multiframe  blind  deconvolution  [38,39]  in  the  polarization  insensitive  case. 
In  the  present  case,  the  joint-likelihood  of  the  combined  polarization  channels  is  pa¬ 
rameterized  by  the  linear  polarization  state  of  the  underlying  scene  and  the  channel 
point  spread  functions.  This  likelihood  function  is  recast  using  the  generalized  ex¬ 
pectation  maximization  (GEM)  algorithm  which,  in  turn,  decouples  the  estimators 
amongst  the  parameters  of  interest. 

Before  proceeding,  it  is  worth  mentioning  some  other  recent  advances  in  po¬ 
larization  imagery  estimation.  In  [51],  a  maximum  a  posteriori  Stokes  parameter 
estimator  is  derived  and  tested.  Unlike  the  present  case,  the  estimator  is  not  in¬ 
tended  to  account  for  sensor  or  atmospheric  point  spread  function  (PSF)  variations 
between  imaging  channels  and,  furthermore,  relies  up  prior  knowledge  that  is  not 
assumed  here.  In  his  dissertation,  Strong  derives  a  blind  maximum  likelihood  esti¬ 
mator  of  a  target  scene  and  a  lumped  polarization  parameter  [43].  This  estimator 
combines  observations  from  one  polarization  channel  and  one  polarization  insensitive 
channel  and  is  loosely  based  on  an  expectation  maximization  approach.  This  differs 
significantly  from  the  present  work  in  that,  here,  the  polarization  state  of  the  target 
scene  is  determined  without  ambiguity  in  degree  or  angle  of  polarization  and  in  that 
any  number  of  polarization  channels  may  be  combined  to  improve  the  estimate. 

The  chapter  begins  with  a  brief  introduction  to  the  GEM  algorithm  and  a 
mathematically  description  of  polarized  light  in  section  4.2.  In  section  4.3,  the 
polarization  image  estimator  is  derived  and  the  case  is  made  for  estimating  channel 
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point  spread  functions  using  Schulz’s  original  method.  For  completeness,  a  brief 
overview  of  the  PSF  estimator  is  also  included.  Finally,  the  algorithm  is  tested  using 
a  laboratory  imaging  polarimeter  in  section  4.4. 

4-2  Preliminary  Notes 

The  purpose  of  this  section  is  to  introduce  and  motivate  the  use  of  the  GEM 
algorithm  for  this  estimation  problem.  In  addition,  a  model  is  produced  for  polarized 
radiation  that  best  compliments  the  GEM  approach. 

4-2.1  The  GEM  Algorithm.  In  maximum  likelihood  (ML)  estimation,  one 
or  more  parameters  of  an  assumed  probability  distribution  for  variable  X  are  esti¬ 
mated  directly  from  observations  of  X.  The  generalized  expectation  maximization 
algorithm  is  an  iterative  method  for  achieving  an  ML  estimate.  The  present  prob¬ 
lem  is  well  served  by  the  GEM  approach  because  an  ML  approach  that  directly 
estimates  the  desired  parameters  is  impractical  to  implement.  The  seminal  paper  on 
the  GEM  algorithm  is  found  at  [1]  and  a  very  readable  overview  at  [29] .  In  what  fol¬ 
lows,  some  basic  steps  are  outlined  to  the  extent  necessary  to  understand  the  GEM 
implementation  in  the  sections  that  follow. 

As  the  name  suggests,  the  GEM  algorithm  consists  of  two  parts:  an  expec¬ 
tation  step  and  a  maximization  step.  In  the  build  up  to  the  expectation  step,  the 
measured  data  (henceforth  referred  to  as  the  incomplete  data)  is  represented  as  the 
aggregate  of  some  underlying  random  variables  (the  complete  data).  The  complete 
data  may  or  may  not  have  physical  significance  but,  in  every  case,  the  aggregated 
complete  data  have  the  same  probability  distribution  as  the  incomplete  data.  In  the 
expectation  step  (E-step),  the  expected  value  of  the  likelihood  of  the  complete  data 
is  determined  when  conditioned  upon  both  the  incomplete  data  and  the  estimate 
from  the  previous  iteration.  In  the  maximization  step  (M-step),  each  estimated  pa¬ 
rameter  is  selected  such  that  it  either  maximizes  the  expected  likelihood  (i.e.  the 
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result  from  the  E-step)  or  at  least  increases  its  value  over  the  previous  iteration. 
The  algorithm  is  stopped  when  the  latest  iteration  of  the  algorithm  meets  some 
predetermined  stopping  criteria. 

In  section  4.3,  Schulz’s  complete  data  for  the  polarization  insensitive  case  is 
decomposed  into  polarized  and  unpolarized  constituents  for  each  polarimeter  imaging 
channel.  As  such,  the  Poisson  distributed  incomplete  data  (photons  counted  at  the 
detector  array)  are  found  to  be  the  sum  of  many  independent  Poisson  distributed 
polarized  and  unpolarized  parts.  In  the  following  section,  the  physical  model  for  this 
decomposition  is  described. 

4-3  Estimator  Derivation 

In  this  section,  the  broad  concepts  introduced  in  section  4.2.1  are  applied  to  the 
polarimetric  blind  deconvolution  problem.  The  section  begins  with  the  incomplete 
data  log  likelihood  to  demonstrate  the  difficulties  encountered  when  attempting  to 
achieve  an  ML  estimate  directly  from  the  probability  distribution  of  the  measured 
data. 


4-3.1  The  Complete  Data  Log-likelihood.  The  incomplete  data  are  the 
collected  images  themselves,  d.  Define  dc(y)  to  be  an  element  in  d  where  the  subscript 
c  specifies  the  channel  in  question  and  y  specifies  a  location  in  the  imaging  space. 
Each  dc(y)  is  modeled  as  an  independent  Poisson  random  variable.  Using  the  same 
convention,  the  mean  of  d  is  i  which  is  a  function  of  the  scene,  o,  and  the  collection 
of  channel  point  spread  functions,  h.  Consequently,  the  mean  of  each  dc(y)  is  ic(y), 
given  by  the  convolution  of  oc  and  hc\ 

ic  (y)  =  °c  (x) hc  (y  ~  (4- !) 

X 

which  is  simply  the  summation  form  of  equation  (2.25).  To  capture  the  variation 
between  channels  due  to  scene  polarization  content,  each  oc  is  decomposed  into  its 
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polarized  constituents: 


X 

oc{x)  =  -Xu{x)  +  Xp(x)  cos2  (a(x)  -  9C)  (4.2) 

Equation  (4.2)  demonstrates  how  polarization  information  is  shared  between 
channels  in  the  log-likelihood.  Now  we  postulate  a  formulation  for  the  polarized 
and  unpolarized  components  of  the  complete  data,  duc(y,  x)  and  dpc(y,  x),  which  will 
alleviate  the  mathematical  difficulties  incurred  by  maximizing  the  likelihood  of  dc(y) 
directly.  To  meet  the  requirements  of  the  GEM  algorithm,  the  complete  data  are 
also  Poisson  random  variables  such  that: 


dc(y)  =  E  due  (y,x)+^dpc(y, 


x) 


(4.3) 


and 


E 


duc(y,  x)  |  =  -A u(x)hc  (y  -  x) 


E 


dpc(y ,  x)  =  \p(x)  cos2  (a(x)  -  9C )  hc  (y  -  x) 


Hence  the  complete  data  log-likelihood,  LCD,  is: 


-A u(x)hc  (y-x) 


(4.4a) 

(4.4b) 


-A u{x)hc  (y-x)  }>  + 


L  (A ui  Ap,  ex,  h )  ^  ^  ^  ^  ^  ^  4  duc(y,  x )  In 

c  y  x  ^ 

\dpc(y,x)  In  [Ap(x)  cos2  (a(x)  —  9C )  hc  ( y  —  x)]  —  \p{x)  cos2  ( a(x )  —  9C )  hc  (y  —  x) 


c  y  x 


(4.5) 


plus  a  term  that  does  not  depend  on  Xu,  Xp,  a  or  hc  and  is  therefore  of  no  consequence 
to  the  estimator. 

4-3.2  The  E-step.  The  expectation  step,  or  E-step,  defines  the  objective 
function  that  is  to  be  maximized  by  careful  selection  of  arguments  in  the  following 
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maximization  step.  Consequently,  there  is  an  obvious  analogy  between  the  role  of 
this  objective  function  and  that  of  the  likelihood  function  in  traditional  maximum 
likelihood  estimation.  For  iteration  n  +  1,  this  objective  function,  Qn+l ,  is  the 
conditional  expectation  of  the  complete  data  log  likelihood.  This  expected  value 
over  the  complete  data  is  conditioned  upon  the  measurements  (i.e.  the  incomplete 
data)  and  the  most  recent  estimates  of  h,  Xu,  \p  and  a: 

Qn+1  (Att,  Xp,  a,h)  =  E  [Lcd  (Xu,  Xp,  a,  h)  \ d,  A",  A”,  a",  hn]  (4.6) 


and,  for  convenience,  define: 


<P£\v,x)=E  dkc(y,x)\dc,\"k,an,K 


(4.7) 


where  k  refers  to  either  u  or  p  depending  on  whether  the  unpolarized  or  polarized 
components  are  under  test.  As  shown  by  [40]  for  a  related  problem,  the  complete 
data  with  these  conditions  are  binomially  distributed  even  though  the  complete  data, 
with  no  conditions  are  Poisson.  The  means  of  these  distributions  are  therefore: 

pc\v ,  x)  =  cos2  ( an(x )  -  6C)  hnc  (y  -  x )  (4.8a) 

^uc\y,x)  =  \^^xni,x)K  (■ y  -  x )  (4-8b) 

The  form  of  Qn+1  may  be  written  out  explicitly  by  substituting  the  (y,  x)  terms 
for  the  dkc  terms  in  (4.5). 

4-3.3  The  M-step.  The  dividends  of  maximizing  Qn+1  in  lieu  of  the  incom¬ 
plete  data  log  likelihood  are  realized  in  this  section.  First,  the  expected  likelihood 
of  the  polarization  components  at  pixel  xq  are  maximized  by  hireling  the  zero  of  the 
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first  partial  derivatives: 


d  Qn+1 
da(xo ) 


dQn+1 

d\p(x0) 


c  y 

d  Qn+1 
d\u(x0) 


^pc\y,x  o) 
Xp(x  o) 


^cos2  (a(x0)  -  9C) 


y-  y-  _  <+ 

c  y  2 


0 


(4.9a) 


(4.9b) 


^pC+1(y,  x0)  tan  (a(x0)  -  9C)  +  ^  Ap(x0)  sin  [2  (a(x0)  -  0C)]  =  0 

c  y  c 

(4.9c) 


where  C  is  the  number  of  available  channels.  Buried  in  each  of  these  results  is  the 
assumption  that  hc(y )  =  1  for  all  c  and  y\  the  conditions  for  this  statement 
will  be  explained  in  section  4.3.4.  At  this  point,  the  proper  selection  of  channel 
orientation  angles,  9Cl  is  critical.  When  6C  is  evenly  distributed  over  all  possible 
linear  polarization  states  (e.g.  0°,  60°,  —60°  for  a  three  channel  system  or  0°,  45°, 
—45°,  90°  in  the  four  channel  case)  the  following  trigonometric  identity  holds  for  all 
possible  a(xo): 

J>s2  (an+1(x0)  -  9C)  =  ^  (4.10a) 

C 

Y,  sin  [2  (an+1(x0)  —  0C)]  =0  (4.10b) 

C 

As  a  result,  equations  (a)  and  (c)  in  (4.9)  are  simplified  and  the  estimators  decouple 
entirely  across  the  parameters.  Hence  the  estimators  for  Xu,  Xp  are  of  the  form: 

Ay+«)  =  §  X  E  w«)  (4.H) 

c  y 

where  k  is  u  or  p  for  the  unpolarized  and  polarized  components  respectively  and  the 
n  +  1  superscript  is  added  to  signify  that  the  estimate  has  been  updated. 

The  reader  may  easily  verify  that  each  7/++1  is  actually  the  fully  polarized 
component  of  image  c.  As  such,  it  is  appropriate  to  introduce  an  intermediate 
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Stokes  vector  Sn+1: 


Sn+1  =  M~HP  (4.12) 

such  that 

=  (4.13) 

where  X1JP  represents  the  vector  of  fully  polarized  components  Equation  (4.13) 

agrees  with  intuition  but  the  proof  is  nontrivial;  the  details  for  the  three  channel 
case  are  worked  ont  in  Appendix  C.  The  second  derivative  test  may  be  applied  to 
the  equations  in  (4.9)  to  verify  that  these  estimates  do  indeed  maximize  Qn+1. 


4-3-4  Incorporating  Schulz’s  GEM  PSF  Estimator.  What  remains  then  is 
to  estimate  the  point  spread  function  for  each  channel.  The  first  partial  derivative 
of  Qn+1  with  respect  to  an  individual  channel  PSF  at  pixel  z  =  y  —  x  is: 


d  Qn+1 
dhc(z ) 


•0p“c+1 


(v,v 


y 


z)  +  4>w1(y,y 

hc(z ) 


*) 


(4.14) 


which  is  identical  to  the  partial  derivative  with  respect  to  the  PSF  of  the  objec¬ 
tive  function  in  the  polarization  insensitive  case.  Consequently,  we  may  directly 
incorporate  Schulz’s  PSF  estimator  (specifically,  in  [39]  section  5)  into  the  present 
case.  What  follows  are  the  salient  points  of  the  original  estimator  with  notes  on 
some  minor  modifications;  the  reader  is  referred  to  the  original  article  for  a  thorough 
derivation  and  an  explanation  of  the  advantages  of  this  approach.  Define: 


h, 


2 


^A(u)exp  [ty>”+1(ii)]  e-,Mux 


(4.15) 


where  A{u)  is  the  aperture  function  of  the  imaging  system  (constrained  such  that 
'YhuA{u)  =  1),  k  is  a  constant  related  incorporating  both  wavelength  and  sampling 
effects,  and  (pP+1  is  the  current  estimate  of  the  atmospherically  induced  phase  func- 


57 


tion  at  the  aperture  for  channel  c.  In  the  case  of  a  circular  aperture  with  radius  r, 
A{u)  =  0  when  ||u||  >  r. 

This  definition  of  hc ,  which  varies  from  Schulz’s  original,  allows  for  ^  hc(y)  = 
1  in  equation  (4.9).  Though  this  is  a  necessary  step  before  successfully  invoking  the 
trigonometric  identities  in  (4.10),  it  comes  with  drawbacks.  In  the  Schulz  model, 
channel  to  channel  gain  variations,  ac,  such  that  E?y  hc(y )  =  ac,  are  accounted  for  in 
the  joint  estimator.  Here,  an  ac  term  would  be  useful  for  modeling  total  transmission 
variations  between  channels.  Total  transmission,  in  this  context,  is  equivalent  to  the 
action  of  an  unknown  neutral  density  filter.  Modification  of  the  estimator  to  allow 
for  these  variations  is  left  to  future  work. 

Since  the  aperture  function  is  known  in  advance,  the  PSF  estimation  problem 
is  recast  into  estimation  of  the  phase  screen,  <pc  via  one  (or  more)  iterations  of  the 
Gerchberg-Saxton  (GS)  phase  retrieval  algorithm  [11]: 


<P 


n+1  _ 

c 


if  E  f  0*0  In  hc{x,  (pc)  >  E  f  0*0  In  hc(x,  <p*), 

X  X 

otherwise 


(4.16) 


where 


f(®)  = 


hc(x,tprc 

Dr 


y 


dc(y) 

(l/)' 


l{y-x) 


given 

Dc  =  Xec+1(t)  =  ^2dc(y) 

x  y 

and  one  iteration  of  the  GS  algorithm: 


(4.17) 


(4.18) 


<Pc  =  ph 


&-1 


VUx^c)  exP  (* '  Ph  (hc(x,  (fc))) 


(4.19) 


where  &  1  is  the  inverse  Fourier  transform  operator  and  “ph”  is  an  operator  that 
extracts  the  phase  angle  from  a  complex  number.  Equation  (4.18)  is  a  statement 


of  energy  conservation.  This  term  does  not  exist  in  its  present  form  in  the  original 
paper;  instead,  there  was  an  additional  constraint  in  Schulz’s  derivation:  ^2xo(x)  = 
1.  The  reader  may  easily  verify,  however,  that  this  change  is  wholly  consistent. 

4-4  Test  Results  from  Laboratory  Data 

In  this  section,  the  polarimetric  blind  deconvolution  algorithm  is  put  to  the 
test  using  data  from  a  laboratory  imaging  polarimeter.  The  test  sensor  consists  of  a 
Photometries  Cascade  512B  camera,  a  single  250  mm  focusing  lens,  and  a  variable 
polarization  analyzer.  The  camera  array  is  512  x  512  pixels  with  a  16  urn  pitch, 
is  cooled  to  —30°  C,  and  has  an  approximately  uniform  response  of  4  photons  per 
count  at  660  nm.  Aside  from  quantization  noise,  the  imager  also  exhibits  “dark” 
noise  and  bias.  At  the  irradiance  levels  shown  below,  this  noise  is  weak  compared  to 
the  photon  dominated  noise  of  the  target.  An  average  dark  bias  is  subtracted  from 
the  data  in  post  processing.  The  lens  is  stopped  down  to  3.175  mm  to  ensure  proper 
sampling.  In  the  configuration  under  test,  the  effective  system  magnification  is  0.22. 

The  test  target  consists  of  two  fully  polarized  parallel  bars,  2  mm  in  length, 
and  back  illuminated  by  a  red  (660  nm  center  wavelength)  diode.  The  polarization 
angles  of  the  two  bars  are  approximately  orthogonal:  2°  for  the  top  bar  and  —83°  for 
the  bottom  bar  (all  angles  are  in  reference  to  the  horizontal  direction  in  the  imagery). 
The  diode  light  passes  through  a  diffuse  screen  prior  to  being  polarized  in  order  to 
even  out  the  illumination  across  the  target. 

The  collected  data  consists  of  three  images,  each  collected  at  a  different  an¬ 
alyzer  orientation:  0°,  60°,  and  —60°.  The  60°  and  —60°  collections  are  corrupted 
by  a  random  plastic  phase  screen  placed  at  the  aperture.  This  plastic  screen  is 
weakly  birefringent.  Between  the  60°  and  —60°  collections,  the  orientation  of  the 
phase  screen  is  rotated.  In  this  way,  each  channel  is  presented  with  a  different  PSF 
analogous  to  the  atmospheric  short  exposure  imaging  case.  Before  processing,  the 
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images  were  cropped  down  to  200  x  200  pixels  and  coarsely  registered.  These  data 
are  shown  in  figure  4.1. 


(a)  the  0°  channel  (b)  the  60°  channel  (c)  the  —60°  channel 

Figure  4.1:  The  estimator  test  data. 


In  the  0°  case,  only  one  bar  appears  in  the  image  because  the  —83°  bar  is 
almost  fully  suppressed  by  the  channel  analyzer.  In  the  other  cases,  phase  error 
dominates.  For  comparison  to  3.1b  and  3.1c,  the  aberration  free  data  for  the  60°, 
and  —60°  channels  are  shown  in  figure  4.2. 


(a)  the  60°  channel  (b)  the  —60°  channel 
Figure  4.2:  The  60°,  and  —60°  channels  in  focus. 


From  the  data  in  figure  4.1  the  algorithm  derives  an  initial  guess  at  the  po¬ 
larized  and  unpolarized  parts  of  the  scene,  as  shown  in  figure  4.3a  and  4.3b.  The 
remaining  images  in  figure  4.3  show  the  results  of  the  algorithm  after  500  iterations 
(at  which  point  the  algorithm  stagnates). 

Recall  that  the  true  target  image  is  fully  polarized.  Consequently,  the  restored 
\u  should  contain  essentially  no  signal;  as  shown  in  figure  4.3c,  this  appears  to  be 
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the  case.  Quantitatively,  Xu  contains  9.34  x  10'  photons  in  the  initial  estimate,  and 
only  1.09  x  105  photons  in  the  final  estimate. 


(a)  initial  Xu  estimate  (b)  initial  Xp  estimate 


(c)  final  Xu  estimate  (d)  final  Xp  estimate 

Figure  4.3:  The  unpolarized  and  polarized  scene  components  before  and  after  restora¬ 
tion. 


Restoration  of  the  fully  polarized  bar  targets  is  achieved.  What  remains  then, 
is  to  consider  the  restored  angle  of  polarization.  From  equation  (4.13),  we  see  that 
an  angle  is  estimated  for  every  pixel  whether  or  not  the  pixel  contains  meaningful 
polarization  content.  As  an  interpretably  aid,  figure  4.4  shows  the  on-target  esti¬ 
mated  angles  masked  off  with  the  image  in  4.3d.  The  estimated  angles  are  —8°  for 
the  top  bar  and  —72°  for  the  bottom  bar. 

The  ultimate  cause  of  this  angle  bias,  which  has  been  observed  across  a  multi¬ 
tude  of  measurements  with  different  phase  screens,  is  elusive.  One  possible  cause  for 
this  discrepancy  is  that  several  real  world  effects  are  not  taken  into  account  in  the 
model.  For  instance,  the  polarizers  are  less  than  ideal  (extinction  ratio  of  ps  106). 
In  addition,  unmitigated  optical  activity,  Fresnel  losses,  or  unmodeled  signal  atten¬ 
uation  at  the  phase  screen  may  be  bias  contributors  as  well.  In  the  latter  case, 
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Figure  4.4:  The  recovered  target  angle  of  polarization  (on  target  pixels  only). 

inclusion  of  a  Schulz  ac  like  term  could  possibly  be  helpful.  It  is  certainly  possible 
that  the  remaining  photons  in  the  unpolarized  image,  if  restored  properly,  could 
possibly  compensate  for  this  bias. 

Finally,  it  is  of  interest  to  present  the  estimated  point  spread  functions  (figure 
4.5).  Sub-figure  4.5a  contains  the  initial  guess  PSF  for  all  channels,  arbitrarily 
selected  to  contain  1  wavelength  of  defocus  aberration.  The  remaining  figures  show 
the  final  estimates  of  the  channel  PSFs  at  the  conclusion  of  the  restoration.  Each 
PSF  in  this  figure  is  magnified  x2  compared  to  the  scale  in  figure  4.1.  Note  that  the 
rotation  of  the  phase  screen  is  evident  between  the  60°  and  —60°  channels.  The  0° 
channel,  which  is  essentially  unaberrated,  tends  toward  a  diffraction  limited  PSF. 

4-5  Technique  Variations:  Prior  Knowledge  and  Multiple  Frames 

Many  variations  of  this  estimator  are  readily  achieved  by  incorporating  priors. 
For  instance,  prior  knowledge  of  the  channel  point  spread  functions  would  be  particu¬ 
larly  useful  in  terms  of  speeding  up  algorithm  convergence  and  increasing  the  overall 
fidelity  of  the  restored  images.  In  this  case,  section  4.3.4  would  be  ignored  and  the 
given  PSFs  could  be  substituted  directly  into  (4.11).  In  a  sense,  this  known  PSF  case 
provides  a  polarimetric,  multichannel  version  of  the  Richardson- Lucy  algorithm. 

There  may  also  be  the  case  where  the  intended  target  is  known  to  be  fully 
polarized  or  unpolarized  (against  a  dark  background).  In  the  fully  unpolarized  case 
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(a)  initial  hc  estimate  (all  (b)  final  0°  PSF  estimate 
channels) 


(c)  final  60°  PSF  estimate  (d)  final  —60°  PSF  esti¬ 
mate 

Figure  4.5:  Close  up  of  the  estimated  point  spread  functions. 

(i.e.  A p{x)  =  0  for  all  x),  the  estimator  reduces  to  Schulz’s  estimator.  In  the  fully 
polarized  case,  the  rule  Xu(x)  =  0  for  all  x  is  enforced  and  only  Xp,  a,  and  h  are 
estimated.  In  more  exotic  variations,  prior  distributions  could  be  imposed  on  the 
angle  or  degree  of  polarization  based  on  source-target-sensor  geometry. 

In  perhaps  the  most  obvious  extension,  multiple  frames  could  be  collected 
for  each  channel.  The  added  noise  and  PSF  diversity  (assuming  the  PSF  is  non¬ 
stationary)  would  “sharpen”  the  likelihood  function  and  could  be  expected  to  im¬ 
prove  the  estimator.  Implementation  of  this  improvement  is  achieved  by  summing 
over  an  additional  index  of  frames  in  (4.11)  with  an  appropriately  modified  per-frame 
weighting  (i.e.  for  N  frames  per  channel,  -AC  instead  of  ^). 

4-6  Improvements  Over  Single  Channel  Deconvolution 

The  results  in  the  previous  section  demonstrate  the  viability  of  the  polarimet- 
ric  blind  deconvolution  algorithm.  What  remains  then,  is  to  determine  how  the 
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algorithm  compares  to  the  alternatives.  Though  this  estimator  has  no  peer  in  the 
literature,  it  is  possible  to  estimate  each  channel  image  and  point  spread  function 
individually  (i.e.  without  knowledge  of  the  other  channels)  and  then  combine  the 
result  in  the  end  to  estimate  Xu,  Xp  and  a  or  any  equivalent  representation.  In  this 
section,  simulated  polarimetric  imagery  is  used  to  compare  the  new  joint  estima¬ 
tor  with  Schulz’s  polarization  insensitive  blind  deconvolution  algorithm  applied  on 
a  per-channel  basis. 

As  a  historical  note,  Schulz  did  not  address  the  problem  of  estimating  an 
image  and  a  point  spread  function  from  a  single  frame  of  data  in  his  original  work. 
However,  in  the  single  channel  case,  the  image  estimator  in  the  Schulz  algorithm 
is  indistinguishable  from  the  Richardson- Lucy  algorithm  which  has  been  applied  to 
the  single  channel  blind  problem  (see,  for  instance,  [10]).  The  difference  between 
these  two  approaches  lies  in  how  the  PSFs  are  estimated.  Since  there  is  no  obvious 
advantage  to  choosing  one  PSF  estimator  over  the  other,  the  Schulz  PSF  estimator 
is  used  here  for  consistency. 

4-6.1  Target  Simulation.  Similar  to  the  laboratory  data,  the  simulated 
target  consists  of  two  fully  polarized  bar  targets  that  are  each  20  pixels  in  length, 
3  pixels  in  width,  and  separated  by  3  pixels.  The  bars  appear  against  a  dark  field 
of  128  x  128  pixels.  Each  illuminated  pixel  contains  104  photons  before  it  is  blurred 
by  the  PSF,  attenuated  by  the  channel  polarizers,  and  corrupted  by  Poisson  noise. 
The  angle  of  polarization  for  each  bar  is  pulled  from  a  uniform  distribution  of  all 
possible  polarization  angles.  The  phase  for  each  channel  point  spread  function  is 
generated  from  a  random  combination  of  the  first  9  Zernike  polynomial  coefficients 
[34],  Each  Zernike  coefficient  is  normally  distributed  with  zero  mean  and  a  variance 
of  1  wavelength  (i.e.  a  phase  of  27r  radians).  The  blurring  and  channel  attenuation 
action  of  the  simulated  sensor  results  in  a  wide  variety  of  image  signal-to-noise  ratios 
even  without  varying  the  illumination  rate.  The  sensor  itself  is  in  the  previously 
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described  0°,  60°,  and  —60°  configuration.  Example  simulated  channels  are  shown 
in  figure  4.6. 


(a)  the  0°  channel  (b)  the  60°  channel  (c)  the  —60°  channel 

Figure  4.6:  Example  simulated  channels. 


4-6.2  Deconvolution  Algorithm  Implementation.  Recall  that  the  polari- 
metric  deconvolution  algorithm  reduces  to  Schulz’s  algorithm  in  the  case  where  the 
target  is  known  to  be  fully  unpolarized  a  priori.  This  is  also  true  when  only  a  single 
channel  image  and  PSF  are  estimated.  In  this  case,  the  channel  images  are  esti¬ 
mated  individually  from  equation  (4.11)  for  Xu  with  c  =  1.  The  PSF  estimator  is  as 
shown  in  section  4.3.4.  After  each  algorithm  iteration  is  complete,  the  Stokes  vectors 
are  formed  from  (2.23).  Similarly,  the  Stokes  parameters  at  each  iteration  can  be 
calculated  from  Xu,  Xpi  and  a  for  each  iteration  of  the  multichannel  estimator: 


S0  = 


Xu  +  Ar 


(4.20a) 


Sr 


Xp 

Vl  +  tan2  a 


{1  -45°  <  a  <  45° 

—1  otherwise 


(4.20b) 


and,  given  this  result,  S2  can  be  calculated  from  (2.13). 

Both  the  single  and  multichannel  algorithms  are  given  the  same  aperture  size 
and  initial  PSF  guess  (arbitrarily  selected  to  be  1  wave  of  defocus).  Both  algorithms 
are  run  for  100  iterations  for  1500  separate  realizations  of  the  target  and  noise. 
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4-6.3  Results.  Normalized  Mean  Squared  Error  ( NMSE )  is  routinely  used 
as  a  image  deconvolution  evaluation  metric  [21].  For  a  true  image  intensity,  /,  and 
an  estimated  image,  /,  NMSE  is  defined  to  be: 


NMSE 


E  f(x)  -  f(x ) 


E/2(^) 


(4.21) 


where  x  is  the  2-dimensional  coordinates  in  the  image  plane.  It  is  clear  from  (4.21) 
that  a  smaller  NMSE  represents  a  better  estimate  of  /.  Typically,  the  NMSE  is 
plotted  against  iteration  number  to  convey  both  the  accuracy  and  rate  of  convergence 
of  the  estimator.  In  the  present  case,  it  is  desirable  to  convey  the  accuracy  and 
rate  of  convergence  for  both  the  image  itself  and  its  polarimetric  properties.  In 
previous  sections,  image  polarization  content  is  conveyed  through  Xu,  \p,  and  a. 
Clearly,  NMSE  applied  to  an  a  “image”  would  be  uninterpretable.  Instead,  the 
Stokes  parameters  (which  all  have  units  of  intensity)  are  fed  into  the  NMSE  equation 
for  consistent  representation.  Finally,  the  goal  is  to  represent  the  aggregate  results 
of  many  simulated  targets;  therefore,  the  following  plots  contain  the  median  (second 
quartile)  NMSE  with  error  bars  representing  the  first  and  third  quartiles  of  these 
data.  Quartiles  are  used  in  lieu  of  mean  and  standard  deviation  to  avoid  error  bars 
with  negative  values  at  low  NMSE.  The  NMSE  is  shown  for  each  Stokes  parameter 
in  figure  4.7. 

It  is  clear  from  these  results  that  the  multichannel  estimator  provides  both 
a  better  median  restoration  of  the  target  (in  all  cases)  and  a  substantially  smaller 
interquartile  range  on  the  restoration  in  the  S\  and  S2  cases.  Recall  that  parameters 
Si  and  S2  carry  all  the  polarization  information  about  the  scene.  Figures  4.7b  and 
4.7c  for  the  single  channel  estimator  show  that,  after  around  25  iterations,  the  NMSE 
actually  increases  slightly  for  S\  and  S'2;  in  other  words,  the  single  channel  algorithm 
produces  a  median  estimate  that  is  getting  farther  from  the  truth  with  additional 
iterations. 
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(a)  So  (b)  Si 


(c)S2 


Figure  4.7:  NMSE  quartiles  for  the  Simulated  Stokes  images. 

Equation  (4.9)  helps  to  explain  the  observed  behavior.  In  the  multichannel 
case,  information  is  shared  between  channels  via  the  estimates  of  Xu,  Xp,  and  a. 
Mathematically,  this  information  transfer  occurs  during  the  sum  over  c  in  (4.9).  The 
result  of  this  transfer  is  a  more  constrained  estimator.  These  image  constraints  also 
indirectly  translate  into  improved  PSF  estimates  even  though  the  PSFs  themselves 
do  not  share  information  across  channels.  In  the  single  channel  case,  without  this 
shared  information,  the  number  of  possible  combinations  of  PSFs  and  images  that 
maximize  Q  increases,  resulting  in  more  estimator  error  on  average.  One  obvious 
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extrapolation  of  this  argument  is  that  the  restoration  would  be  improved  further  if 
additional  channels  were  added  to  the  estimator. 

f.7  Further  Analysis  of  the  Angle  Bias  Issue 

In  this  section,  simulation  is  used  to  determine  whether  or  not  the  linear  polar¬ 
ization  angle  error  shown  in  section  4.4  is  a  symptom  of  an  actual  bias  in  a  statistical 
sense.  The  simulated  data  is  as  described  in  section  4.6.  In  this  case,  the  multichan¬ 
nel  deconvolution  algorithm  is  run  for  200  iteration  over  1500  realizations  of  the 
data.  (200  iterations  are  used  is  this  case  to  ensure  near  convergence  of  the  results.) 
The  restored  target  polarization  angles,  a,  are  taken  to  be  the  average  over  the  true 
bar  target  area  and  the  bias  is  defined  as: 

b  =  a  —  a  (4.22) 

where  b  is  defined  between  —90°  and  90°  (example:  a  —  a  =  120°  =>  b  =  —60°).  The 
simulation  results  show  an  average  bias  of  0.5°  with  a  standard  deviation  of  16.3°. 
Consequently,  even  in  the  absence  of  calibration  and  unmodeled  noise  errors,  the 
laboratory  results  in  section  4.4  are  typical,  at  least  in  the  sense  that  the  reported 
angle  errors  are  within  a  standard  deviation  of  the  mean  error  in  this  simulated  case. 
A  histogram  of  the  angle  bias  results  (bin  width:  3.66°)  is  shown  in  figure  4.8. 


Figure  4.8:  Results  of  the  angle  bias  simulation. 


4-8  Chapter  Summary 

In  this  chapter,  a  maximum  likelihood  multichannel  blind  deconvolution  esti¬ 
mator  is  developed  from  the  generalized  expectation  maximization  algorithm.  The 
efficacy  of  the  algorithm  is  demonstrated  using  test  data  from  a  three  channel  imaging 
polarimeter  in  a  laboratory  setting  and  in  simulation.  The  laboratory  reconstructed 
test  imagery  is  very  close  to  the  true  underlying  target  scene  in  terms  of  the  es¬ 
timated  total  polarized  and  unpolarized  content  though  somewhat  biased  in  terms 
of  angle  of  polarization.  Though  it  is  difficult  to  identify  the  cause  precisely,  this 
type  of  error  has  been  observed  across  several  such  experiments  with  different  phase 
screens.  In  simulation,  the  polarization  angle  estimator  is  shown  to  be  unbiased 
under  ideal  image  acquisition  conditions.  Possible  error  sources  include  unmodeled 
noise  and  calibration  errors  in  the  target,  phase  screen,  or  sensor. 
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V.  Detection  through  Obscurations  Using  Optimized  Temporal 

Polarization  Imaging 

In  this  chapter,  temporal  fluctuations  in  polarization  signature  are  maximized  by 
forming  linear  combinations  of  polarization- sensitive  intensity  channels.  These 
combination  channels,  which  are  the  result  of  a  principal  components  analysis,  are 
also  reduced  in  rank  and  therefore  optimized  for  further  processing.  Temporal  fluc¬ 
tuations  are  emphasized  as  a  means  to  detect  polarizing  objects  under  substantial 
but  weakly  polarizing  obscurations  that  would  otherwise  preclude  detection  in  a  tra¬ 
ditional  intensity  image.  The  theoretical  work,  which  is  an  expansion  of  previously 
published  research,  is  demonstrated  for  three  and  four  channel  polarization  imaging 
systems.  Finally,  the  theory  is  tested  through  simulated  examples  using  an  empirical 
target  obscuration  model. 

5.1  Target  Detection  in  Obscurations 

Several  remote  sensing  techniques  have  been  applied  to  the  problem  of  target 
detection  in  the  presence  of  heavy  obscurations.  Recent  examples  of  modes  applied  to 
this  problem  include  radar  [30],  hyperspectral  imaging  [17],  and  L1DAR  [37]  among 
others.  The  present  work  takes  the  first  steps  toward  addressing  the  obscured  target 
detection  problem  by  exploiting  temporal  fluctuations  in  target  polarization  signa¬ 
ture  due  to  time  varying  obscurations.  What  follows  is  a  description  of  a  methodol¬ 
ogy  believed  to  optimize  a  polarization  imaging  system  to  tackle  this  obscured  target 
problem. 

Consider  a  passive  incoherent  imaging  system  with  any  number  of  distinct 
channels.  Each  channel  is  most  responsive  to  a  different  orientation  of  linear  polar¬ 
ization.  Under  these  constraints,  a  sensing  modality  is  developed  to  detect  polarizing 
targets  through  weakly  polarized  or  unpolarized  obscurants.  A  realistic  example  of 
this  scenario  would  be  detection  of  vehicles  under  heavy  forest  canopy  from  an  air- 
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borne  polarization  imaging  system.  Canopy  motion  or  the  movement  of  the  aircraft 
will  provide  the  sensor  with  small,  transient,  polarized  reflections  off  the  target. 
Taken  individually  (i.e.  from  a  single  frame),  these  transient  events  may  be  statisti¬ 
cally  indiscernible  from  the  obscurants  in  the  space  around  them.  Now  consider  an 
ensemble  of  these  frames,  all  registered  to  the  same  scene.  The  obscured  target  pix¬ 
els  should  have  an  appreciably  larger  variance  than  those  pixels  filled  with  obscured 
natural  background  (which  is  also  largely  unpolarized).  The  premise  of  this  research 
is  realized  by  determining  the  configuration  of  the  polarization  imager  that  is  best 
suited  to  detect  this  target  variance  in  polarization  state  while  suppressing  other 
sources  of  scene  fluctuation  (for  instance,  variance  in  total  intensity)  and  correlation 
between  channels. 

J.S.  Tyo  published  a  method  for  optimizing  a  multi-channel  polarization  imag¬ 
ing  system  using  principal  components  analysis  [45] .  His  work  provides  the  founda¬ 
tion  for  this  research,  but  not  without  modification.  Specifically,  Tyo’s  derivation  re¬ 
quires  all  incident  radiation  to  be  monochromatic.  In  doing  so,  this  work  is  applicable 
to  all  permutations  of  clliptically  and  linearly  polarized  light  but  not,  unfortunately, 
to  unpolarized  light,  which  is  of  specific  interest  here  (unpolarized  light  is,  among 
other  things,  polychromatic  [42]).  To  reconcile  this  deficiency,  section  5.2  applies  this 
Stokes  model  to  Tyo’s  optimization  strategy.  This  result  is  then  extended  further 
to  define  the  optimum  combination  of  channels  for  detecting  temporal  fluctuations 
in  polarization  state.  Section  5.3  applies  this  strategy  to  two  notional  polarization 
imaging  systems  including  one  of  Tyo’s  examples,  which  is  shown  to  be  equivalent 
in  this  new  formulation.  Section  5.4  describes  an  empirical  obscured  target  model 
which  is  employed  in  a  simulated  case  study. 

5.2  Channel  Optimization 

As  in  each  of  the  previous  chapters,  the  polarization  imaging  system  is  com¬ 
posed  of  N  channels,  each  of  which  is  most  responsive  to  a  different  orientation  of 
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linearly  polarized  light.  Optimized  combinations  of  these  channels  are  generated 
(either  through  hardware  or  software  implementation)  and  passed  onto  an  image 
accumulator  where  a  to-be-determined  number  of  frames  are  stored  until  statisti¬ 
cal  relevance  is  achieved.  In  this  case,  the  stored  images  are  then  used  to  form  a 
standard  deviation  image  for  each  of  the  optimized  combination  channels. 

Equipped  with  a  Stokes  description  of  intensity  at  the  detector  array,  Tyo’s 
method,  which  is  presented  below,  can  be  employed  to  construct  a  set  of  channel 
combinations  that  are  uncorrelated,  possibly  reduced  in  rank,  and  hence  simplified 
for  further  processing.  These  uncorrelated  combination  channels  are  then  further 
optimized  for  maximum  response  to  temporal  fluctuations  in  polarization  state. 


5.2.1  Tyo’s  Method.  The  total  intensity  measured  for  each  channel  is  rep¬ 
resented  by  an  N  x  1  vector,  I.  The  correlation  matrix  for  an  N  channel  polarization 
imaging  system  is  given  by: 

C  =  E^[I-IT]  (5.1) 


where  E^[-]  is  the  expected  value  taken  over  all  possible  polarization  states. 
Following  Tyo’s  prescription,  decorrelation  of  these  channels  is  achieved  by  solving 
for  a  matrix  X  such  that 


XTCX 


A  k  i  =  j 
0  i^j 


(5.2) 


where  A&  are  the  N  eigenvalues  of  C  and  i ,  j  are  the  row  and  column  coordinates  of 
XTCX.  In  other  words,  XTCX  is  a  diagonal  matrix  whose  diagonal  entries  are  the 
eigenvalues  of  C .  In  this  arrangement,  these  eigenvalues  are  also  the  second  moment 
of  the  (yet  to  be  formed)  uncorrelated  combination  channels.  This  fact  will  be  useful 
later  on.  To  achieve  the  required  solution  for  (5.2),  the  columns  of  X  are  taken  to 
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be  the  N  normalized  eigenvectors  of  C,  Vjp 

A"  =  [vi...vN]  (5.3) 

Having  determined  X,  the  transformed  channels  (eigenchannels) ,  Z,  are  given 
by: 

Z  =  XTI  (5.4) 

The  channels  Z  are  uncorrelated  over  the  ensemble  of  possible  linear  polarization 
states  thus  an  expansion  of  Tyo’s  original  work  is  achieved.  What  remains  is  to  find 
a  subset  of  Z  that  is  most  sensitive  to  temporal  fluctuation  in  polarization  state. 

5.2.2  Temporal  Polarization  Expansions.  If  this  were  a  true  principal 
components  analysis,  the  next  step  would  be  to  select  a  subset  of  Z  to  form  a  basis 
of  possibly  reduced  rank  to  represent  all  of  I.  Selection  for  this  new  basis  would  be 
determined  from  the  magnitude  of  the  associated  eigenvalues  [28].  Instead,  channels 
in  Z  are  selected  based  on  their  ability  to  maximize  signal  variance  due  to  temporal 
fluctuations  in  degree  of  polarization.  More  succinctly,  if  7,  some  minimum  allowable 
variance  threshold,  is  set,  then  an  eigenchannel  in  Z  is  retained  if: 

E$  EP  ( Zi  -  Zi)2  >7  (5.5) 

where 

Zi  =  EP  [Z^  (5.6) 

where  the  subscript  P  indicates  expectation  over  degree  of  polarization.  Recalling 
that 

A.  =  E(  [Z?]  (5.7) 

equation  (5.5)  can  be  simplified  by  rearranging  the  order  of  expectations: 
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Ep  [Aj]  —  E g 


(5.8) 


>  7 


This  expectation  rearrangement  is  acceptable  because  there  is  no  reason  to  believe 
that  angle  of  polarization  and  degree  of  polarization  are  statistically  dependent  quan¬ 
tities. 


Methods  for  choosing  7  can  vary,  especially  if  a  priori  knowledge  of  the  dis¬ 
tribution  in  P  is  available.  Without  this  knowledge  of  P  however,  it  will  be  shown 
that  selecting  7  =  0  is  sufficient  in  several  useful  cases. 


5.3  Example  Polarization  Imaging  Sensors 

The  previously  described  methodology  is  completely  general.  The  following 
examples  contain  a  number  of  simplifications  that  are  not  required  to  gain  a  solu¬ 
tion  but  do  allow  for  more  compact  results.  Each  polarizer  is  assumed  to  be  ideal. 
Additionally,  it  will  be  assumed  that  all  incident  radiation  is  linearly  polarized, 
unpolarized  or  partially  linearly  polarized.  In  other  words,  there  is  no  elliptically 
polarized  light  and  S', 3  =  0.  Next,  it  is  assumed  that  the  distribution  in  total  intensity 
is  independent  from  the  distribution  of  linear  polarization  states  (this  assumption  is 
always  true  for  unobscured  targets). 

Since  there  is  now  only  linear  polarization  to  consider,  the  angle  of  linear 
polarization,  £,  will  be  defined  as  in  Chapter  II.  Finally,  it  is  assumed  that  the 
polarized  channels  are  evenly  spaced  on  the  same  interval  as  £.  In  addition  to  be 
being  a  simplifying  step,  this  assumption  also  ensures  that  one  of  the  eigenchannels 
in  Z  is  the  sum  of  all  the  original  intensity  channels  (as  shown  by  Tyo). 

5.3.1  The  3- Channel  Case.  The  3-channel  case  is  presented  here  for  two 
reasons;  the  first  is  because  it  reduces  the  number  of  steps  required  to  attain  infor¬ 
mative  results  and,  second,  it  provides  an  example  of  how  the  derivation  in  [45]  is 
equivalent  to  the  case  presented  here. 
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Intensity  detected  by  our  3-channel  imaging  polarimeter  is  given  by  substitu¬ 
tion  into  equation  (2.8): 


Io  =  y  (1  +  Pcos2£) 


T  _  u 

-*60  —  ^ 


T  _  u 

-*-60  —  — 


cos2£  —  \^3  sin  2£ 


cos2£  +  y/3  sin  2£ 


(5.9a) 

(5.9b) 

(5.9c) 


where  the  subscripts  0,  60  and  —60  indicate  the  angle  of  polarization  (in  degrees) 
preferred  for  transmission  through  the  channel  polarizer.  To  reiterate  a  previous 
point,  the  orientations  of  these  filters  are  not  arbitrary,  they  are  selected  such  that 
they  are  evenly  spaced  on  the  interval  of  £. 


Given  i  and  j  as  indices  of  the  row  and  column  entries  in  the  3x3  correlation 
matrix: 


71-/2 

Cij  I  Iiljdi ; 


7 r 


(5.10) 


— 7r/2 


such  that  i  or  j  —  1  corresponds  to  Jo,  i  or  j  =  2  corresponds  to  Iqq,  and  so  on.  The 
matrix  C  has  eigenvalues 


and  (normalized)  eigenvectors 


A  -  A  -  3p23° 

-  A3  -  -^g- 


(5.11) 


1 

_ 1 

I 

l 

1 

’7  =  * 

1 

v*=7i 

0 

^3  =  76 

-2 

1 

1 

1 

(5.12) 
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Note  that  several  other  eigenvectors  with  eigenvalues  equal  to  A2  are  possible,  but 
all  are  equivalent.  This  eigenvector  result  is  in  agreement  with  the  special  case  of 
the  3-channel  derivation  in  [45]. 


Moving  on,  the  eigenchannels  are  given  by  substitution  of  (5.9)  into  (5.3)  and 


(5.4): 


Z  = 


^75  (h  +  -^60  +  J-6Cl) 

7S  (J-60  -  Jo) 

77g  (Jo  —  2/60  +  I- 60 ) 


(5.13) 


Since  there  is  no  a  priori  knowledge  about  how  these  eigenchannels  are  distributed 
with  respect  to  P,  the  test  from  equation  (5.8)  will  be  executed  with  7  =  0.  First, 
note  that 


Zi  —  Ep 


77^  (Jo  +  ho  +  J-60) 


=  E„ 


V3So 

2 


=  V3S0 
2 


(5.14) 


hence,  for  the  first  uncorrelated  channel  at  least,  the  exact  distribution  of  the  degree 
of  polarization  is  irrelevant.  This  result  can  be  readily  substituted  into  the  test  from 
equation  (5.8): 

"  Vas0'  2~ 


3So 


Ef 


=  0 


(5.15) 


Interestingly,  even  though  this  channel  would  be  the  largest  contributing  principal 
component  in  a  traditional  sense  (by  virtue  of  its  dominating  eigenvalue),  it  fails  to 
exceed  the  7  =  0  variance  threshold  and  therefore  must  be  discarded. 


The  eigenvalues  of  the  two  remaining  channels  are  equal  and  both  dependent 
on  P  so,  for  any  meaningful  distribution  on  P,  their  variance  will  exceed  a  threshold 
of  7  =  0.  Consequently,  both  of  these  channels  are  retained  as  part  of  the  scheme  to 
maximize  temporal  variation  in  degree  of  polarization. 


5.3.2  The  4-Channel  Case.  A  4-channel  imaging  polarimeter  (with  equally 
spaced  polarization  channels)  is  of  interest  because  three  of  the  Stokes  parameters 
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can  be  calculated  directly  using  addition  and  subtraction  on  combinations  of  the 
channel  outputs  [5]: 

•So  =  4  +  4o 

Si  =  Iq  —  4o  (5.16) 

S-2  =  I45  —  1 135 

In  this  case,  the  eigenvalues  of  the  correlation  matrix  are  given  by: 


(5.17) 


(5.18) 


Now  the  variance  threshold  test  in  (5.8)  is  applied  to  Z.  No  additional  manipulation 
is  required  to  show  that  the  sum  channel  (Ai  =  Sq)  will  not  pass  the  threshold 
test  for  the  same  reasons  that  the  sum  channel  did  not  pass  in  the  3-channel  case. 
Furthermore,  the  channel  corresponding  to  A2  =  0  contains  no  useful  information 
and  would  be  discarded  by  both  a  traditional  principal  components  analysis  and  the 
temporal  variance  maximization  analysis  presented  here. 

The  channels  corresponding  to  A3  and  A4  remain.  Observation  of  (5.18)  shows 
that  these  channels  are  actually  just  scaled  versions  of  the  Stokes  parameters  S2  and 
,Sj  as  seen  in  (5.16).  This  result  may  suggest  that  it  is  only  required  to  cast  the 
measurements  into  S2  and  Si  to  maximize  sensitivity  to  fluctuations  in  degree  of 
polarization.  Proof  of  this  assertion  is  left  to  future  work. 
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5-4  Analysis  Using  Simulated  Data 

The  purpose  of  this  section  is  to  demonstrate  the  utility  of  the  imaging  scheme 
derived  in  section  5.3  using  simulation.  The  simulation  consists  of  a  two  plane 
obscured  target  model  that  varies  in  time  in  a  statistically  predictable  way.  The  test 
data  are  presented  along  with  the  results  from  each  eigenchannel. 


5-4-1  Obscured  Target  Model.  The  simulation  consists  of  a  simple  empir¬ 
ical  obscuration  model  for  a  Stokes  vector  image  in  the  object  plane  and  an  image 
plane  model  found  using  an  arbitrary  point  spread  function.  Each  point  in  the  object 
plane,  S'(u,  n ),  can  be  decomposed  into  two  basis  planes:  f(u,  n ),  the  front  or  obscu¬ 
ration  plane,  and  b(u,  n),  the  back  or  target  plane.  Here,  u  is  the  two-dimensional 
coordinates  in  the  object  plane  and  n  is  the  frame  number.  Each  point  in  both  /  and 
b  is  described  by  a  4-dimensional  Stokes  vector  corresponding  to  the  case  where  each 
basis  plane  is  directly  illuminated  by  the  source.  Whenever  the  scene  is  illuminated 
by  natural  light,  each  point  in  S'  can  be  represented  as  a  linear  combination  of  / 
and  b.  Specifically: 


S'(u,  n )  =  t(u,  n)f(u,  n)+ 
[1  —  t(u,  n )]  r(u,  n)b(u,  n) 


(5.19) 


where  /  and  b  are  scaled  by  t(u,  n ),  the  fraction  of  a  pixel  obscured  by  /,  and  r(u,  n), 
the  fraction  of  the  maximum  possible  b  that  is  actually  reflected.  At  each  position  u, 
S'  is  a  realization  of  a  random  process  defined  by  the  four  underlying  processes:  /, 
b ,  t  and  r.  Processes  /  and  b  are  each  assumed  to  have  two  components.  The  larger 
component  varies  slowly  with  the  respect  location  of  the  illumination  source  (i.e. 
the  sun)  while  the  smaller  component  varies  rapidly  with  small  changes  in  surface 
orientation.  The  slowly  varying  components  of  /  and  b  are  treated  as  constants. 
Conversely,  t  and  r  vary  with  each  frame  and  depend  partially  on  the  density  of  the 
obscurant.  In  cases  where  the  obscurant  is  a  natural  media  (e.g.  a  forest  canopy), 


78 


the  obscurant  is  assumed  to  be  weakly  polarized  when  compared  to  the  underlying 
target. 

Having  previously  established  a  requirement  for  incoherent  illumination,  the 
Stokes  image  at  the  detector  is  given  by 

S(x,  n )  =  [h®  5']  (x,  n)  (5.20) 

where  h  is  the  sensor  specific  point  spread  function.  Equation  (5.20)  can  be  readily 
derived  for  each  entry  in  S.  Consider  the  formation  of  an  incoherent  image,  /, 
without  regard  to  the  polarization  of  the  incident  light: 

I  =  h®S'0  (5.21) 

where,  borrowing  from  the  definition  of  the  Stokes  vector,  S'0  is  the  total  intensity 
of  the  scene  as  predicted  by  geometrical  optics.  S'0  can  be  decomposed  into  parts 
that  represent  orthogonal  polarization  states  and,  since  convolution  is  distributive, 
so  can  /: 

I  =  h®(l±  +  -T||)  =  h  ®  /_l  +  h  <g>  7||  (5.22) 

By  careful  selection  of  the  orthogonal  components  in  this  equation  (refer  to  equation 
(5.16)),  one  can  form  any  of  the  other  terms  in  the  Stokes  vector.  For  instance: 

Si  =  h  <g>  Iq  —  h  <g)  /go  —  h<S>  S[  (5.23) 

Similarly,  the  same  process  can  be  used  to  define  S2  in  terms  of  S'2  and  so  on  for  S3. 

Equation  (5.19)  is  the  simplest  possible  model  for  mixing  two  images  and 
does  not  include  polarization  effects.  Instead,  polarization  effects  are  added  in  later 
through  judicious  selection  of  Stokes  vectors  for  different  classes  in  the  scene. 
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Additionally,  the  model  in  (5.19)  contains  no  noise  terms.  Though  not  neces¬ 
sary  to  demonstrate  the  concepts  from  the  previous  section,  a  realistic  noise  model 
would  be  critical  for  establishing  limits  on  the  effectiveness  of  this  obscuration  de¬ 
tection  scheme. 

Furthermore,  equation  (5.19)  ignores  the  effects  of  thickness  and/or  depth  in 
the  obscuration  “plane”.  To  illustrate  the  limitations  on  this  model,  consider  the 
case  where  the  primary  obscuration  is  much  nearer  to  the  sensor  than  the  target. 
In  this  case,  the  obscurations  would  effect  the  formation  of  the  average  h(x)  as  well 
as  S'.  These  effects  on  h{x)  are  well  described  in  [12].  At  a  minimum,  the  ratio  of 
the  target- to-obscuration  distance  over  sensor-to-target  distance  should  be  <<  1  for 
model  applicability  in  this  sense.  More  specific  bounds  on  model  applicability  are 
left  to  future  work. 

5-4-2  Test  Data  Sets.  As  described  in  the  previous  section,  the  obscured 
target  model  requires  a  front  plane  (figure  5.1),  containing  the  primary  obscura¬ 
tion,  and  a  back  plane  figure  (5.2),  which  contains  the  targets  and  some  natural 
background. 


Figure  5.1:  The  front  plane  or  foliage  image  used  during  the  simulation. 

For  the  purposes  of  simulation,  these  images  are  broken  down  into  three  pri¬ 
mary  classes:  targets,  foliage  and  grass.  Obviously,  more  classes  exist  in  both  of  these 
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Figure  5.2:  The  unobscured  back  plane  image  containing  the  targets  used  during  the 
simulation. 

images.  These  additional  classes  have  been  ignored  for  expediency  and,  regardless, 
the  ideas  put  forward  in  section  5.3  can  be  tested  without  them. 

Each  class  is  assigned  an  average  Stokes  vector: 


*S 'grass 


*S foliage 


Si 


target 


1  -.05  .05  0 
1  .03  .03  0 
1  -.46  0  0 


T 

T 

T 


(5.24) 


Each  of  the  grass  and  foliage  classes  have  an  average  degree  of  polarization  of  only  a 
few  percent  while  the  target,  which  is  made  to  represent  a  near  specular  reflection, 
is  significantly  higher. 


The  average  percent  pixel  obscuration,  E[t(u,n)],  is  99.9%  with  a  standard 
deviation  of  0.03%.  In  other  words,  even  though  the  front  plane  is  “flat”  it  almost 
completely  obscures  the  back  plane  at  all  times.  Likewise,  the  average  percentage  of 
the  maximum  possible  reflection  off  the  backplane,  E[r(u,  n)]  is  50%  with  a  standard 
deviation  of  10%.  The  purpose  of  r(u,n )  is  to  partially  capture  the  effects  of  depth 
in  the  obscuration  by  forcing  the  backplane  to  be  only  partially  illuminated  with 
significant  overall  fluctuation. 
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The  point  spread  function  used  to  blur  the  composite  Stokes  image  is  shown  in 
figure  5.3  in  the  same  scale  as  the  image.  These  parts  come  together  under  equation 
(5.20)  to  form  simulated  Stokes  images  in  the  detector  plane. 


Figure  5.3:  Simulated  PSF. 


5. 4-3  Results.  All  simulated  results  are  shown  with  a  100%  linear  stretch 
and  the  MATLAB  default  color  map.  Figure  5.4  shows  a  single  simulated  image 
from  the  Z\  through  Z±  channels  given  by  equation  (5.18).  Though  present  in  the 
data,  the  locations  of  the  obscured  targets  are  imperceptible. 

True  to  the  predicted  eigenvalues  for  the  4-channel  case  in  equation  (5.17);  Z\ 
appears  to  contain  the  most  information  (largest  eigenvalue),  Z-2  contains  virtually 
no  information  (0  eigenvalue),  while  Z3  and  Z4  both  contain  some  information  (but 
still  no  indication  of  the  targets).  Equation  (5.24)  shows  that  the  targets  are  most 
highly  polarized  in  ,S'i  so,  if  the  targets  could  be  detected,  it  would  likely  occur  in 
eigenchannel  Z4  (figure  5.4d).  Regardless  of  channel,  target  detection  through  the 
provided  obscuration  using  only  a  single  frame  appears  unlikely. 

Figure  5.5  contains  the  standard  deviation  images  of  channels  Z\  through  Z4. 
Each  standard  deviation  image  is  calculated  from  25  separate  realizations  of  an 
eigenchannel  image.  Each  additional  image  is  generated  with  a  different  randomized 
t(u,n )  and  r(u,n).  Perfect  registration  between  images  is  assumed  throughout. 
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Figure  5.4:  A  single  frame  of  the  simulated  eigenchannels. 


The  derivation  in  section  5.3  for  the  4-channel  case  predicts  that  channels  Z\ 
and  Z2  should  be  discarded  because  their  variance  is  0  with  respect  to  P,  which  is 
assumed  to  the  variable  driving  fluctuation  from  image  to  image.  Figures  5.5a  and 
5.5b,  the  standard  deviation  images  of  Z\  and  Z2,  appear  to  generally  uphold  this 
prediction  with  some  interesting  exceptions. 

First,  note  the  foreground  grass  in  figure  5.5a.  No  portion  of  the  background 
scene  penetrates  this  portion  of  the  image  and,  consequently,  the  variance  across  this 
region  is  very  low  (the  regional  variance  is  low  but  not  zero  because  of  point  spread 
function  effects).  In  other  words,  since  the  grass  here  is  impenetrable,  it  violates 
the  original  assumptions  made  about  the  obscuration  and  hence  its  behavior  in  the 
standard  deviation  image  is  not  predictable  by  this  model. 

Second,  there  is  some  indication  of  the  right-most  target  in  the  Z\  standard 
deviation  image  even  though  this  optimization  method  predicts  that  there  should  be 
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(c)  z3  (d)  Z4 

Figure  5.5:  Standard  deviation  images  of  the  eigenchannels. 


none.  To  explain  this  apparent  incongruity,  recall  that  the  distribution  over  all  pos¬ 
sible  intensities  was  assumed  to  be  separable  from  the  distribution  of  possible  linear 
polarization  states.  This  assumption  is  required  by  equation  (5.10)  for  formation  of 
the  multi-channel  covariance  matrix  but  is  violated  later  when  the  obscured  target 
model  is  defined.  Additionally,  variance  that  occurs  purely  due  to  fluctuation  in  total 
intensity  is  ignored  throughout  this  development.  Though  clearly  this  assumption 
is  not  entirely  accurate,  it  still  appears  to  be  a  reasonable  approximation,  especially 
when  the  information  content  in  this  channel  is  compared  with  that  in  the  Z3  and 
Z4  channels. 

For  the  same  reasons,  the  Z 2  channel  and  corresponding  standard  deviation 
image  also  contain  some  information  content  though  none  of  it  is  apparently  useful. 

Little  further  commentary  is  required  to  show  that,  as  predicted  by  this  opti¬ 
mization  method  with  7  =  0,  channels  Z3  and  Z4  clearly  contain  the  most  informa- 
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tion  about  the  obscured  targets.  Consequently,  no  target  information  would  have 
been  lost  if  channels  Z\  and  Z<±  had  been  discarded  completely. 

Up  to  this  point,  no  consideration  has  been  given  to  the  original  (i.e.  correlated) 
intensity  channels.  Figures  5.6  and  5.7  contain  a  single  frame  and  standard  deviation 
image  for  the  J90  channel,  which,  for  the  targets  as  defined,  should  contain  the  most 
target  signal. 


Figure  5.6:  A  single  simulated  frame  from  the  /go  channel. 


Figure  5.7:  The  J90  standard  deviation  image. 


The  /90  standard  deviation  image  does  show  some  evidence  of  the  targets 
though  the  quality  of  the  information  present  is  far  inferior  to  that  contained  in 
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the  Z4  standard  deviation  image.  For  brevity’s  sake,  the  other  intensity  channels  are 
not  shown  because  the  results  are  much  the  same  as  the  /90  channel. 

5. 5  Summary 

A  principal  components  analysis  of  a  polarimetric  imaging  system  is  derived 
for  the  case  of  polarized,  partially  polarized,  and  unpolarized  radiation.  This  result 
is  then  applied  to  the  problem  of  imaging  through  obscurations  by  discarding  the 
eigenchannels  that  are  likely  to  be  insensitive  to  temporal  fluctuations  in  polariza¬ 
tion  state.  The  optimized  temporal  polarization  imaging  method  presented  in  this 
research  provides  a  viable  scheme  for  target  detection  in  the  presence  of  obscura¬ 
tions.  Defense  of  this  assertion  is  provided  in  the  theoretical  development  and  tested 
against  a  simulated  dataset. 

In  the  simulated  example,  with  no  prior  knowledge  of  the  target,  this  method 
reduces  the  number  of  channels  required  for  effective  target  detection  from  4  to  2 
with  relatively  trivial  pre-  and  post-processing  requirements.  This  notional  sensor 
was  configured  such  that  three  of  the  Stokes  parameters  can  be  calculated  directly 
from  the  measured  intensities.  Under  these  circumstances,  the  channel  reduction 
result  is  expected  since  only  two  of  the  three  Stokes  parameters  depend  on  degree 
of  polarization.  Referring  back  to  the  derivation  at  (5.18),  perhaps  the  more  inter¬ 
esting  result  from  this  example  is  that  Stokes  parameters  are  also  shown  to  be  the 
optimum  combination  channels.  Finally,  this  example  is  a  special  case;  the  demon¬ 
strated  method  allows  for  calculation  of  the  optimum  combination  channels  for  any 
polarization  imaging  sensor. 

While  the  results  are  compelling,  the  simulation  ignores  three  important  diffi¬ 
culties  that  will  arise  in  real  world  applications.  First,  image  to  image  registration 
is  assumed  to  be  perfect  for  all  channels  at  all  times.  Second,  real  targets  in  this 
scenario  may  be  obscured  in  part  by  objects  that  are  completely  opaque  at  all  times, 
effectively  breaking  up  the  target  in  the  standard  deviation  image.  Third,  all  signals 


were  assumed  to  be  noiseless.  Compensation  for  these  difficulties,  and  more  minor 
unnamed  others,  will  be  fleshed  out  in  future  work  using  held  data. 
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VI.  Conclusion 


The  three  preceding  chapters  address  problems  involving  polarimetric  image 
registration,  restoration,  and  analysis  in  the  presence  of  noise  and  a  non- 
deterministic  channel.  Unlike  the  bulk  of  the  related  literature,  these  problems 
are  approached  from  an  estimation  theory  perspective  combined  with  a  complete 
physical  model  of  polarization.  As  is  demonstrated,  this  combination  is  fertile  with 
research  opportunities  that  draw  upon  research  in  traditional,  polarization  insensi¬ 
tive  image  processing.  In  this  section,  the  previous  chapters  are  summarized  and 
concluding  remarks  are  presented. 

In  chapter  III,  the  Cramer-Rao  lower  bound  on  image  registration  errors  by 
Robinson  and  Milanfar  is  generalized  to  the  case  of  multichannel  polarimetric  im¬ 
agery.  In  doing  so,  the  problem  is  posed  as  a  joint  estimation  of  both  the  translational 
registration  errors  and  of  the  polarimetric  image  itself.  The  primary  difficulty  in  cal¬ 
culating  the  bound  directly  is  that  the  result  requires  inversion  of  an  enormous  (even 
by  the  standards  of  desktop  computers)  Fisher  information  matrix.  This  inversion 
is  made  tractable  by  applying  matrix  theory  to  the  FIM.  In  doing  so,  the  bounds 
largely  decompose  into  components  that  are  easily  identifiable  in  the  processing 
chain:  target,  receiver,  and  estimator. 

The  bound  is  also  used  to  describe  three  and  four  channel  polarimetric  imaging 
systems  and  the  special  case  of  N  polarization  insensitive  images.  For  N  frames  of 
polarization  insensitive  imagery,  the  Cramer-Rao  bound  derived  by  Robinson  and 
Milanfar  is  shown  to  be  the  limit  of  the  bound  as  N  becomes  large.  In  the  three 
channel  polarization  case,  the  Cramer-Rao  bound  for  any  joint  estimator  is  infinity. 
Consequently,  the  registration  parameters  and  Stokes  images  must  be  estimated 
separately.  An  additional  treatment  of  joint  estimation  versus  external  measurement 
can  be  found  in  appendix  B.3.  Finally,  the  form  of  the  bound  suggests  that  the 


optimum  channel  spacing  is  60°  and  45°  respectively  in  the  three  and  four  channel 
cases. 

In  chapter  IV,  an  iterative  maximum  likelihood  blind  deconvolution  algorithm 
is  derived  using  expectation  maximization.  This  algorithm  leverages  a  complete 
data  model  based  on  the  polarized  and  unpolarized  components  of  the  scene  and  a 
uniform  spacing  between  channels  to  decouple  the  joint  estimator  into  four  separate 
estimators  of  the  polarized  intensity,  angle  of  polarization,  unpolarized  intensity, 
and  point  spread  function.  The  resulting  algorithm  is  shown  to  be  viable  using 
imagery  collected  in  the  laboratory.  Simulated  data  is  used  to  demonstrate  the 
improvement  gained  using  this  multichannel  approach  when  compared  to  the  only 
viable  alternative,  blind  deconvolution  of  each  channel  individually.  Finally,  the 
estimated  errors  in  polarization  angle  from  the  laboratory  data  are  shown  to  be 
typical  of  results  found  in  simulation. 

In  chapter  V,  a  principal  components  analysis  of  polarization  imagery  is  com¬ 
bined  with  temporal  variation  thresholding  for  polarimetric  imagery.  The  goal  of 
this  effort  is  to  reduce  the  number  of  channels  required  to  detect  potential  targets 
in  the  presence  of  obscurations.  Results  are  demonstrated  theoretically  for  the  three 
and  four  channel  cases.  In  the  three  channel  case,  the  PCA  optimized  channels 
are  shown  to  agree  with  previous  research  in  which  only  fully  polarized  radiation 
is  considered.  In  the  four  channel  case,  the  PCA  and  temporal  variation  optimized 
channels  are  shown  to  be  the  Stokes  parameters  S\  and  S2 .  This  four  channel  case 
is  then  validated  using  simulation  and  a  simple  random  obscuration  model. 

6.1  Research  Extensions 

The  Cramer-Rao  bound  research  opens  the  possibility  for  substantial  future 
work.  This  research  could  follow  the  lead  of  Robinson  and  Milanfar  [33]  and  address 
the  problem  of  registration  estimator  bias  for  a  specific  algorithm  as  it  applies  to 
polarimetric  imagery.  Additionally,  Yetik  and  Nehorai’s  work  could  be  expanded  to 


include  a  polarimetric  bound  for  other  types  of  image  deformation  (e.g.  rotation, 
shearing,  etc.)  or  bounding  estimator  error  for  polarimetric  registration  algorithms 
that  use  feature  extraction  techniques.  Additionally,  a  search  could  be  conduced  for 
a  joint  estimator  that  meets  the  Cramer- Rao  bound.  Finally,  the  idea  of  representing 
bounds  on  image  registration  as  a  realization  of  underlying  parameters  (in  this  case, 
as  the  weighted  sum  of  the  Stokes  parameters)  may  be  used  to  derive  bounds  for 
registering  multi-modal  images  if  a  similar  underlying  relationship  exists. 

In  addition  to  the  estimator  variations  presented  in  section  4.5,  the  blind  de- 
convolution  research  should  be  expanded  to  include  sensor  noise  and  the  effects  of 
non-ideal  polarization  analyzers.  The  angle  bias  question  may  possibly  be  addressed 
by  systematically  eliminating  calibration  issues  or  incorporating  them  into  the  es¬ 
timator.  Compensation  for  between  channel  (neutral  density)  transmission  differ¬ 
ences  is  also  desirable.  Additionally,  the  algorithm  is  readily  expanded  to  included 
sensitivity  to  elliptical  states  of  polarization.  Finally,  it  may  be  of  interest  to  the 
astrophysical  community  to  reevaluate  polarization  imagery  from  the  Hubble  Faint 
Object  Camera  in  the  context  of  this  new  tool. 

The  optimized  analysis  research  can  be  expanded  in  a  number  of  directions. 
First,  the  PCA-thresholding  research  must  be  tested  against  field  data  to  account 
for  real-world  variations  in  targets  and  obscurations.  Such  a  test  would  be  relatively 
easy  to  conduct  and  far  superior  to  simulation.  The  reasons  for  this  necessity,  if  not 
obvious,  are  discussed  in  section  5.5.  Second,  the  distribution  of  possible  polarization 
states  (used  to  calculate  the  channel  correlation  matrix)  can  be  improved  to  include 
knowledge  of  the  sensor  location  and  possibly  the  source  with  respect  to  the  target 
area.  In  one  possible  scenario,  the  location  of  the  aircraft  mounted  sensor  and 
the  position  of  the  sun  can  be  used  to  approximate  likely  polarization  states  for 
“flat”  dielectric  or  metal  targets  underneath  trees.  The  physical  mechanisms  for 
this  distribution  are  discussed  qualitatively  in  chapter  II.  Finally,  an  actual  target 
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detector  (in  the  classical  sense)  can  be  added  to  the  output  of  the  optimized  channels 
to  possibly  eliminate  or  reduce  the  need  for  a  human  analyst. 
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Appendix  A.  Useful  Background  Derivations 

The  purpose  of  this  appendix  is  to  expand  on  some  of  the  intermediate  results  from 
Chapter  II. 


A.l  The  Field  Along  the  Transmission  Axis  of  an  Ideal  Retarder- Polarizer  Pair 

At  first  glance,  the  result  in  equation  (2.2),  which  is  the  scalar  sum  of  orthog¬ 
onal  electric  field  components  after  passing  through  a  retarder  and  polarizer,  may 
appear  to  be  in  error.  Here,  the  result  is  shown  to  be  accurate  using  the  Jones 
matrices  for  polarized  light.  A  full  accounting  of  the  Jones  calculus  can  be  found  in 
many  general  optics  texts  (for  instance:  [3,16,31]).  Jones  matrices  describe  direct 
manipulation  of  the  electric  field  and  are  therefore  quite  useful  for  building  up  to 
the  definition  of  the  Stokes  parameters.  After  this  point,  the  Jones  matrices  are 
abandoned  and  the  Stokes  parameters  are  manipulated  directly. 

The  vector  representation  of  the  field  in  (2.1)  is  given  by: 


ux  (t ) 
uy  (t ) 


(A.l) 


At  the  sensor,  this  field  passes  through  a  retarder  (or  waveplate),  which  imposes 
a  phase  shift  of  f  between  the  x  and  y  components,  and  then  a  polarizer,  which 
has  a  primary  transmission  axis  along  angle  9.  Each  of  these  optical  elements  is 
represented  by  a  2  x  2  Jones  matrix: 
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The  resulting  field,  E',  can  be  expressed  in  a  rotated  coordinate  system  (using  an¬ 
other  Jones  matrix)  where  one  of  the  principal  axes  is  also  the  polarizer  transmission 
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axis: 

cos  9  sin  9 
—  sin  6  cos  8 

- N/ - 

rotation 

such  that,  in  this  new  coordinate  system: 

ux  (t)  e ^  cos  9  +  uy 

0 

which  is  the  result  in  equation  (2.2). 

A. 2  The  Decomposition  of  Partially  Linearly  Polarized  Light 

The  decomposition  of  partially  linearly  polarized  light  into  its  constituents 
is  of  great  importance  in  Chapter  IV.  Consequently,  the  result  in  equation  (2.15) 
is  derived  here  in  detail.  Partially  polarized  light  with  Stokes  vector,  S,  can  be 
represented  as  the  sum  of  two  Stokes  vectors:  Su,  the  unpolarized  component,  and 
Sp,  the  fully  polarized  component. 

S  =  Su  +  Sp  (A. 5) 


where 

1 
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where  ip  is  the  angle  of  polarization  (equation  2.13)  and  P  is  the  degree  of  polarization 
in  (2.12).  From  equation  (2.7),  the  measured  intensity  for  the  unpolarized  component 
is  therefore: 


/«  =  (!-  P)S0  =  K  (A. 7) 

and,  likewise,  the  measured  intensity  for  the  polarized  component  is: 

Ip  =  -PS0  (1  A  cos 2'0 cos 29  +  sin 2^ sin 29) 

=  l,PSo  [1  +  cos2(^  -  9)}  (A.8) 

=  PS0  cos2(-0  —  9) 

and  so,  by  defining  \p  =  PSo,  the  total  intensity  at  the  detector  is  given  by: 

I  —  lu  +  Ip  —  2^u  A  cos2  —  P)  (A. 9) 

which  is  equation  (2.15). 
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Appendix  B.  More  on  the  Polarimetric  Cramer- Rao  Bound 

The  first  two  sections  of  this  appendix  are  devoted  to  developing  mathematically 
necessities  for  the  bound  derivation  in  Chapter  Ill.  The  remaining  three  sections  are 
devoted  to  expansions  on  this  theme  that  were  not  part  of  the  originally  published 
work. 


B.  1  More  on  the  V  and  V  Matrices 

The  matrices  V  and  V  play  a  central  role  in  the  bound  calculation.  In  addi¬ 
tion,  the  V  matrix  provides  the  connection  between  the  previous  work  by  Robinson 
and  Milanfar  [33]  and  the  generalization  provided  here.  Consequently,  some  useful 
properties  of  these  matrices  are  developed  in  this  section. 

First,  there  is  an  important  connection  between  the  V,  V,  and  : 


(V  \T  v  _  ajiajk 

{- n-ij )  nkj  -  ~r~ 
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and 


(Hij)T  Hki  = 
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for  all  i  A  3- 


Second,  each  v*  is  established  in  the  plane  of  the  image.  Therefore,  a  common 
direction  vector,  x,  for  all  relevant  derivatives  can  be  defined  via  the  chain  rule: 


d_ 


(B.3) 


Consequently,  each  entry  in  V  and  V  is  defined  in  a  common  coordinate  system. 
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B.2  Detailed  Derivation  of  the  Stokes  Parameter  Bound 

The  average  bound  can  be  determined  by  calculating  the  trace  of  equation 
(3.21)  via  application  of  the  following  properties  of  the  trace  and  the  Kronecker 
product.  For  the  trace  [13]: 


tr  (A  +  G)  =  tr  (A)  +  tr  (G)  (B.4a) 

tr  ( CGD )  =  tr  ( DCG )  (B.4b) 

tr  ( EF )  =  vec  ( ET ) T  vec  (F)  (B.4c) 

where  A  and  G  are  square  matrices  and  C,  D,  E,  and  F  are  any  matrices  such  that 
CGD  and  EF  are  square  matrices.  Operator  tr  represents  the  trace  and,  for  any 
matrix  A,  vec(A)  is  defined  to  be  an  ordered  stack  of  the  columns  of  A.  For  the 
Kronecker  product  (also  from  [13]): 

vec  ( AGC )  =  ( CT  ®  A)  vec  (G)  (B.5a) 

(A  ®  C)  (G  ®  D)  =  {AG  ®  CD)  (B.5b) 

both  of  which  always  hold  whenever  AGC,  AG,  and  CD  are  defined.  Finally,  note 
that,  in  contrast  to  regular  matrix  multiplication: 

{A  0  G)t  =  At  ®Gt  (B.6a) 

{A  •  G)t  =  At.Gt  (B.6b) 
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Armed  with  the  preceding  formulae,  derivation  of  the  average  bound  for  each  Stokes 
image  is  straightforward.  Starting  with  the  second  term  in  equation  (3.32): 


tr  [kb"1  ( HBvHt )  =  tr  (HBvHt)\ 

=  vec  (<S>f^)Tvec  (. HBVHT ) 

=  vec  ($“T$-1)T  (H  ®  H )  vec  (. Bv ) 

=  vec  [Ht  H]T  vec  (Bv) 

Fortunately,  expression  HT  H  is  identical  in  form  to  equation  (3.22)  and 

can  therefore  be  simplified  in  the  same  manner: 


Ht  (4>rT4>r‘)  H  =  <7- Ht  (C^C-'  ®  Ip,,p>)  H 

=  a^Wsi  •  (v  +  v) 

where 

WSi  =  M  (C'-TC'-1)  Mt  ®  l2x2 
Invoking  the  implicit  symmetry  of  Wgi,  V,  and  V : 

tr  ( HBvHt )  $~T]  =  a2 vec  [wSi  •  (v  +  y)] T  vec  (Bv) 

=  a2tr  [wSi  •  (v  +  y) 

This  result  can  be  folded  back  into  (3.31)  to  produce: 

m  =  p  [tr  (r-1)]  +  pr  .  (V  +  V)  Bv 

2 

=  a2c^  +  °^tr  •  (y  +  y) 

which  is  identical  to  equation  (3.35). 


(B.8) 


(B.9) 


(B.10) 


(B.  11) 


97 


B.3  External  Measurement  Versus  Joint  Estimation 


The  bounds  derived  in  III  provide  the  minimum  achievable  variance  for  an 
unbiased  estimator  in  each  of  the  following  cases: 

•  estimates  of  translational  misregistration  when  the  underlying  Stokes  parame¬ 
ters  are  known 

•  estimates  of  the  Stokes  images  when  the  registration  parameters  between  chan¬ 
nels  are  known 

•  joint  estimates  of  registration  and  Stokes  parameters 

What  remains  then  is  to  specify  a  bound  for  the  cases  where  misregistration  and 
Stokes  estimates  are  determined  separately.  This  bound  can  be  compared  to  the 
joint  estimation  bound  to  determine  which  of  these  paths  yields  the  best  estimator. 

For  the  unrelated  problem  of  bounding  estimates  on  dispersive  wave  param¬ 
eters,  Kimball  (et  al.)  [19]  provides  a  definition  for  the  total  error  achieved  when 
Cramer- Rao  bounds  are  considered  in  conjunction  with  errors  in  externally  provided 
parameters.  Like  in  the  present  work,  Kimball  employs  the  natural  partitioning  of 
the  Cramer-Rao  bound  using  block  matrix  methods.  Going  a  step  further,  Kimball 
then  replaces  the  Cramer-Rao  bound  for  a  subset  of  the  parameters  with  a  externally 
measured  covariance  matrix.  In  this  way,  errors  in  externally  provided  parameters 
can  be  propagated  into  the  expression  for  errors  in  the  estimates  of  the  remaining 
parameters. 

To  apply  Kimball’s  ideas  to  the  present  problem,  define  e2  to  be  an  externally 
determined  covariance  matrix  for  the  channel  registration  parameters.  Note  that  e2 
may  be  less  than  V  if  the  misregistration  between  images  has  a  strong  deterministic 
component.  The  “total  error”  matrix  for  the  Stokes  image  estimates  is  determined 
by  replacing  V  with  e2  in  equation  (3.21): 


4  =  S"1  +  S~lHelHTS~l 


(B.12) 


which,  in  turn,  may  be  simplified  using  the  matrix  methods  presented  in  section 
3.4.3: 

(4i>  =  ^c-1  +  4«r  [n-s. .  (V  +  C)  4]  (B.13) 

where  the  i  subscript  indicates  which  Stokes  parameter  is  under  test. 

Equation  (B.13)  provides  the  means  for  evaluating  the  impact  of  alternate 
forms  of  misregistration  correction  on  Stokes  image  estimation.  This  analysis  can 
be  used  when  considering  design  trade-offs  in  conjunction  with  other  factors  such  as 
cost,  processing  power,  weight,  and  so  on. 

B.4  Interpreting  the  CRLB  Using  Correlations 

Assume  that  the  collected  images  are  bandlimited,  sampled  at  the  Nyquist 
frequency  or  higher,  and  periodic.  Under  these  conditions  Robinson  and  Milanfar 
show  that,  in  the  two  image  case,  the  terms  in  the  matrix  V  are  independent  of 
the  translational  error  between  the  images.  Proof  of  this  fact  is  accomplished  in 
the  spatial  frequency  domain  via  Parseval’s  theorem.  This  result  can  be  readily 
generalized  to  the  present  N  channel  case.  Consider  Parseval’s  theorem  applied  to 
a  non-zero  submatrix  in  V : 


—  OO 


(B-14) 

where  Ft(uj)  is  the  Fourier  transform  of  image  /)  at  frequency  coordinates  u.  The  jeo 
terms  in  the  right  hand  side  are  a  result  of  differentiation  in  the  left  hand  side  and 
the  terms  are  the  Fourier  shift  theorem  representation  of  the  registration 

errors.  Unlike  the  two  image  case,  it  is  not  immediately  clear  from  this  result  that 
the  Fisher  information  matrix  is  independent  of  the  unknown  translations  between 
the  images,  and  v/.  To  resolve  the  issue,  the  integral  on  the  right  hand  side  of 
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(B.14)  can  be  rearranged, 


1 

2-ko2 


J  LULUTFk  (u)  F* 


—  OO 


(B.  15) 


and  treated  as  an  inverse  Fourier  transform  evaluated  at  the  origin  of  x: 


2ll(72 


L0C0TFk  (u)  F*  (ca)e-nvfe~v*  )ue~jx  Uduj 


(B.16) 


x=02xl 


Bringing  this  all  together: 


1 


a 


2 


-ld_ 

a2  <9x 


[fk  *  fi)  (x) 


T 


x=vfc-v; 


(B.  17) 


where  *  is  the  correlation  operation.  Evaluation  of  the  correlated  images  at  x  = 
Vk  —  vi  is  equivalent  to  evaluation  of  the  registered  images  at  the  origin;  in  other 
words,  the  Fisher  Information  matrix  does  not  depended  on  the  translation  between 
images  and  the  generalization  of  the  two  image  case  is  confirmed. 

The  observation  that  V  and  V  can  be  calculated  from  the  second  order  deriva¬ 
tives  of  the  cross-correlation  is  also  of  interest.  First,  correlation  is  not  a  one-to-one 
mapping  over  all  possible  image  pairs.  As  such,  two  unique  ensembles  of  images  have 
the  same  Fisher  information  matrix  if  their  cross-correlations  are  the  same.  Indeed, 
the  cross-correlations  need  only  to  be  equal  near  x  =  Vk  —  vj  since  only  the  second 
derivative  behavior  at  this  point  is  of  consequence.  Heuristically,  this  result  means 
that  images  with  similar  second  order  statistics,  for  instance  imagery  over  similar 
terrain,  will  result  is  similar  bounds. 
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For  the  joint  Stokes  restoration  problem,  the  cross-correlation  between  images 


can  be  recast  in  terms  of  the  Stokes  vectors: 


-5  S 


EE  0  /;.  i  0  l  j  (  S  j  'k  Sj 


(B.18) 


'  '  \  L  i  ]  J  / 

It  is  therefore  clear  that  polarization  channel  orientation  may  serve  to  suppress  or 
reenforce  the  contribution  to  the  correlation  amongst  the  individual  Stokes  images. 


B.5  A  CRLB  for  the  Blind  Deconvolution  Polarization  Parameterization 

As  demonstrated  throughout  this  dissertation,  it  is  often  useful  to  switch  be¬ 
tween  parameterizations  of  polarization  state  in  order  to  best  suit  a  particular  prob¬ 
lem.  In  Chapter  IV,  the  Stokes  vector  parameterization  was  replaced  with  an  in¬ 
tensity  and  angle  parameterization  in  order  to  meet  the  requirements  of  the  com¬ 
plete  data  model.  As  a  bridge  between  these  two  chapters,  this  section  describes 
a  mechanism  for  transforming  the  Stokes  parameterized  CRLB  into  one  using  the 
intensity-angle  parameterization. 

Recall  that  the  original  Stokes  image  parameter  vector  has  the  form: 


9  = 


(B.19) 


where  vectors  vx  are  registration  parameters  and  Sx  are  the  true  Stokes  parameter 
values  for  each  point  in  the  image.  Similarly,  a  likelihood  parameterization  must  be 
defined  for  the  intensity-angle  parameterization: 


w  = 


V  .  .  .  V 


N 


\T  \T 

Au  Ap 


a 


(B.20) 


where  the  vx  vectors  are  the  same  as  before  and  each  of  Xu,  Xp,  and  aT  represent 
the  unpolarized  intensities,  polarized  intensities,  and  angles  of  polarization  for  each 
point  in  the  image.  Note  that  both  6  and  w  have  the  same  dimensions.  Using 
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these  parameter  vectors,  the  transformed  Fisher  information  matrix,  ./',  is  defined 
by  Scharf  [36]  as: 

./'  =  GJGt  (B.21) 

where  J  is  the  FIM  in  equation  (3.1)  and  the  elements  of  G  are  given  by: 

=  H  <B22> 

where  i ,  j  are  the  row  and  column  indices  of  G.  As  before,  the  CRLB  in  the  unbiased 
case  is  given  by  the  inverse  of  J'. 

What  remains  then  is  to  define  the  entries  in  G.  First,  note  that  the  registration 
terms  are  the  same  in  both  parameterizations;  in  that  case: 


dvi 


{hx2  if  i  =  j, 
02x2  if  i  7 ^  3 


(B.23) 


which,  together  over  i,j  —  2  to  N ,  form  a  submatrix  of  G,  henceforth  referred  to  as 
9i- 

Next  we  note  that  the  translation  parameters  and  the  intensity  parameters 
(in  either  form)  are  not  dependent  on  each  other.  Thus  the  partial  derivatives  of 
any  of  these  parameters  with  respect  to  any  translational  parameter  is  zero  (e.g. 
din  =  xp2)  =  O2 xP2i  etc)- 

The  remainder  of  the  terms  describe  the  transformation  between  the  Stokes 
and  intensity- angle  parameters.  At  each  point  mra,  the  Stokes  parameters  can  be 
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described  as: 


S0(mn)  =  Xu(mn)  +  Ap(mn)  (B.24a) 

=  ====T(o:(mn))  (B.24b) 

a/1  +  tan  a( mn) 

S2 (mra)  =  S'i(mn)  tan(2o:(mn))  (B.24c) 


where 


T(a) 


1  -45°  <  a(mn)  <  45° 

—  1  otherwise 


(B.25) 


which,  all  together,  can  be  used  to  populate  the  remaining  components  of  G.  A 
number  of  these  derivative  terms  can  be  expressed  analytically.  Assuming  n  =  k ,  for 
the  unpolarized  component: 


d»So(mn)  _  ^ 

<9Au(mfc) 

dSi( mw)  _  dS2(nin)  _  Q 
d\u(mk)  dXu(mk) 


(B.26a) 

(B.26b) 


the  polarized  component: 


cASVifm., 


=  1 


dXp(mk) 
dSi(mn)  T  [a(m„)] 


(9Ap(mfc)  yA  +  tan2  a(mn) 
dS2{  mn)  aS'i(mn) 

aV^  =  ^W«tan[2“(m")1 


(B.27a) 

(B.27b) 

(B.27c) 


and  the  polarization  angle: 


dSp  (m  w)  _  Q 
da(mk) 


(B.28) 


Note  that  and  d^2j'mn}  do  not  have  a  simple,  closed  form  solution.  Finally,  in 

oaym^)  oa(m 17  j  i 

all  cases  with  n  /  k,  the  partials  in  B.26,  B.27,  and  B.28  are  all  zero.  This  region  of 
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G  is  collectively  referred  to  as  r/2.  Thus: 


G  = 


9 1 

03p2x2(iV-l) 


02(Af-l)x3p2 

92 


(B.29) 
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Appendix  C.  Derivation  of  the  polarization  angle  estimator  for  the 

three  channel  case 

Here  we  derive  the  polarization  angle  estimator  for  the  case  of  polarization  analyzers 
oriented  at  0°,  60°,  and  —60°.  Three  trigonometric  identities  are  required: 


.  „  .  tan  a  —  tan  9C 

tan  (a  -  6C)  =  - - - 

1  +  tan  a  tan  tq 

(C.la) 

tan  2a  =  ( - hi)  tan  a 

\  cos  2a  ) 

(C.lb) 

2  tan  a 

sin  2a  =  - 2 - 

tan  a  +  1 

(C.lc) 

Apply  the  first  identity  to  equation  (4.9)  part  (c),  divide  out  the  common  denomi¬ 
nator,  and  evaluate  the  result  at  each  9C: 


v^3  (4/ 3  —  \l/2)  [tan2  a  +  l]  +4  (dq  +  T2  +  d/3)  tana  —  3dq  tana  [tan2  a  +  l]  =0 

(C.2) 

where  for  compactness,  dq  =  Jfiyrfpk1(yixo)  and  ol  =  an+1(a;o).  This  result  is 
reduced  further  by  applying  (4.12): 


S%+1  =  l  (*  i  +  ^2  +  *3)  (C.3a) 

S?+1  =  |  (d>i  -  d>2  -  d>3)  (C.3b) 

S^1  =  (*2  -  *3)  (C.3c) 


such  that  (C.2)  becomes: 


-?S9n+1  +  35iq+1  2tana  +  5  (SZ+1  +  S[l+1)  tana  =  0. 


tan2  a  +  1 


(C.4) 
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Now  apply  (2.14),  along  with  the  remaining  identities  to  achieve: 

q  q  crn+1 

--S?+1  +  3Sn+1  +  tan  2a  =  0 
2  2  2  2 

which,  when  solved  for  a,  becomes  equation  (4.13). 


(C.5) 
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Appendix  D.  Practical  Laboratory  Polarimetric  Imagers 

Laboratory  collected  imagery  plays  an  important  role  in  the  bound  examples  from 
chapter  III  and  in  the  validation  of  the  deconvolution  algorithm  in  Chapter  IV. 
For  completeness  and  as  an  aid  to  future  progression  of  this  work,  the  laboratory 
system  are  described  in  detail  in  this  appendix.  The  constructed  sensors  fall  into 
two  categories  based  on  how  the  data  are  collected:  serial  (each  channel  is  collected 
sequentially)  and  parallel  (all  channels  are  collected  simultaneously).  The  tradeoffs 
between  these  configurations  will  be  discussed  in  the  sections  that  follow.  Obviously, 
there  are  many  practical  polarimetric  imaging  collection  schemes  other  than  those 
addressed  here;  these  particular  configurations  represent  two  simple  and  effective 
approaches  that  can  be  obtained  with  minimal  cost  and  time  investment. 

D.l  The  Serial  Imager 

The  serial  imager  consists  of  a  variable  analyzer,  focusing  lens,  aperture  stop, 
and  focal  plane  array  (field  stop).  The  basic  setup  is  shown  in  figure  D.l.  The  re¬ 
quired  sequence  of  images  is  collected  by  manipulating  the  analyzer  between  frames. 


Figure  D.l:  Diagram  of  the  single  channel  polarimetric  imager. 

The  variable  polarization  analyzer,  aperture  stop,  and  focusing  lens  are  all  co¬ 
located.  The  front  and  back  of  this  assembly  are  shown  in  figure  D.2.  The  variable 
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analyzer  determines  the  orientation  of  the  preferred  angle  of  polarization  with  respect 
to  some  reference  (in  this  case,  the  plane  of  optical  bench).  The  preferred  orientation 
is  adjusted  by  mechanically  rotating  the  polarizing  element  between  images. 


(a)  the  target  side  (b)  the  camera  side 

Figure  D.2:  The  serial  imager  analyzer-stop- lens  assembly. 


For  a  convex  focusing  lens  with  focal  length  /,  the  thin  lens  equation  gives  us: 


111 
7  —  _  +  v 
J  o  i 


(D.l) 


where  o  is  the  distance  from  the  focusing  lens  to  the  target  and  /  is  the  focal  length 
of  the  lens.  Overall  target  magnification,  rriT  is  given  by: 


i 


rriT  =  ~ 
o 


(D.2) 


In  this  scenario,  the  image  at  the  focal  plane  will  be  real  and  inverted. 

The  aperture  stop  is  located  between  the  variable  analyzer  and  the  focusing 
lens.  In  this  scenario,  the  primary  purpose  of  the  aperture  stop  is  to  ensure  proper 
sampling  of  the  point  spread  function.  Somewhat  counter-intuitively,  proper  sam¬ 
pling  of  the  PSF  in  the  image  plane  is  ensured  by  only  considering  the  sampling  in 
the  aperture.  In  the  quasi-monochromatic  Fresnel  diffraction  regime,  sampling  in 
the  aperture  plane,  A u,  is  related  to  sampling  in  the  frequency  space  of  the  image, 
A£,  and  is  given  by  [9]: 
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A  u  =  A(Xz 


X  z 

NAx 


(D.3) 


where  A  is  the  center  wavelength  of  the  illumination,  z  is  the  distance  between  the 
exit  pupil  and  the  imaging  plane,  N  is  the  number  of  samples  across  the  imaging 
plane,  and  Ax  is  the  sampling  in  the  image  plane  (fixed  by  choice  of  the  FPA).  It  is 
worth  noting  that,  unless  the  optics  in  question  are  very  large  or  z  is  very  short,  this 
imager  will  always  be  operated  in  the  Fresnel  regime.  Furthermore,  it  is  assumed 
that  the  FPA  is  square.  In  this  setup,  the  exit  pupil  is  always  the  aperture  stop  itself 
and  z  —  i  in  figure  D.l. 

For  incoherent  image  formation,  the  PSF  is  proportional  to  the  magnitude 
squared  of  the  Fourier  transform  of  the  aperture  function  (see  equation  4.15).  From 
[9],  the  maximum  spatial  frequency  passed  by  this  optical  system  (i.e.  the  cutoff 
frequency)  in  this  scenario  is: 


C.=  £  (D.4) 

where  R  is  the  radius  of  the  aperture  stop.  Hence,  to  exceed  the  Nyquist  threshold 
for  sampling,  NA(  >  2(c  or,  in  other  words: 


R  < 


X  z 
AAx 


It  is  also  often  convenient  to  express  this  result  as  N  >4^ 


(D.5) 


Finally,  the  aperture  stop  and  focusing  lens  together  also  determine  the  irra- 
diance  that  reaches  the  FPA  (throughput).  According  to  [31],  the  irradiance  at  the 
FPA  is  proportional  to  j (i.e.  the  inverse  of  the  f-number  squared).  For  a  fixed 
focal  length,  reducing  the  aperture  radius  by  a  factor  of  2  means  that  only  |  of  the 
photons  reach  the  FPA.  Therefore,  the  camera  integration  time  must  be  4  times  as 
long  to  achieve  the  same  photon  (and  noise)  levels. 
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Table  D.l:  A  workable  serial  polarimetric  imager. 


/ 

0 

i 

R 

Ax 

25  cm 

29  cm 

133  cm 

1.588  mm 

16  um 

The  choice  of  focusing  lens  and  aperture  stop  is  a  trade-off  between  the  space 
available  for  the  laboratory  setup,  the  desired  target  magnification,  the  sampling 
requirements  imposed  by  the  FPA,  and  throughput.  The  polarization  analyzer  itself 
plays  no  particular  role  in  this  analysis.  As  a  guideline,  the  greatest  system  flexibility 
is  achieved  by  choosing  the  largest  possible  aperture  given  the  available  FPA,  lenses 
and  lab  space.  Magnification  is  a  secondary  issue  (unless  the  image  exceeds  the  space 
on  the  FPA)  and  one  may  always  reduce  the  size  of  the  aperture  without  violating  the 
sampling  requirements  in  (D.5).  An  additional  benefit  to  selecting  larger  apertures 
is  that  larger  apertures  are  easier  to  aberrate  (which  is  useful  primarily  for  testing 
blind  deconvolution  algorithms).  Table  D.l  contains  the  design  used  to  produce  the 
laboratory  example  in  Chapter  IV. 

D.1.1  Advantages  and  Disadvantages.  Besides  ease  of  setup,  the  other 
advantages  of  the  serial  imager  include  precise  and  easily  verifiable  angular  channel 
positioning  without  speciality  machined  parts.  Additionally,  each  channel  may  be 
focused  and/or  stopped  individually  with  minimal  extra  effort.  Assuming  the  same 
focal  plane  is  used  in  both  configurations,  more  pixels  are  available  per  channel 
which  translates  into  improved  resolution  or  increased  field  of  view  depending  on 
the  specific  configuration. 

There  are  also  disadvantages  to  the  serial  imager.  First,  the  variable  polarizer 
must  be  adjusted  between  collections  for  each  channel.  This  substantially  increases 
time  spent  in  the  laboratory  and  often  requires  that  collected  images  be  registered 
in  post  processing.  Second,  the  target  must  remain  static  between  collections. 
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D.2  The  Parallel  Imager 

The  imaging  polarimeter  described  in  this  section  captures  both  target  in¬ 
tensity  and  polarization  state  simultaneously.  The  optical  system  consists  of  three 
primary  components:  the  channel  lens  array,  an  image  reducing  optic,  and  field  stop. 
The  aperture  stop  and  channel  analyzers  are  co-located  with  the  lens  array.  These 
components  are  arranged  as  shown  in  figure  D.3. 


Field  Stop 
(at  the  target) 


Reducing 

Optics 


Focusing  Lens  Array, 
Analyzers, 
and  Aperture  Stop 
(co-located) 


Figure  D.3:  Diagram  of  the  parallel  polarimetric  imager. 


The  target  scene,  defined  by  the  limits  of  the  field  stop,  is  imaged  separately 
by  each  of  the  lenses  in  the  channel  lens  array.  Each  channel  is  most  sensitive  to 
polarization  in  a  different  way.  Consequently,  the  intermediate  images  formed  by  the 
lens  array  communicate  different  target  polarization  information.  The  intermediate 
images  are  spatially  separated  from  each  other  by  the  diameter  of  the  lenses  them¬ 
selves.  This  composite  of  images  is  made  to  fit  onto  the  focal  plane  array  via  the 
reducing  optic.  Once  digitized,  the  composite  image  is  processed  to  form  a  single 
polarization  image. 

The  channel  lens  array  is  an  optical  component  that  produces  three  separate 
images  of  the  target  scene  each  of  which  responds  to  polarized  light  in  a  different 
way.  Each  of  three  channels  consist  of  an  aperture  stop,  polarizing  filter  (analyzer) 
and  a  small  diameter  lens.  Baffling  is  placed  around  the  array  and  in  the  small 
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region  in  the  array  center  to  reject  unintended  illumination  from  the  target.  Each 
channel  is  identical  except  for  the  orientation  of  the  analyzer  (0°,  60°,  and  —60°). 
The  radius  of  the  stop  is  determined  again  by  (D.5)  to  meet  sampling  requirements, 
though  calculation  of  z  is  somewhat  more  complicated  (as  shown  below).  The  front 
and  back  of  this  array  are  shown  in  figure  D.4. 


(a)  the  target  side  (b)  the  camera  side 

Figure  D.4:  The  parallel  imager  channel  lens  assembly. 


The  three  images  formed  by  the  lens  array  alone  (along  with  the  empty  space 
between  images)  are  much  too  large  to  fit  on  a  typical  CCD  array  therefore  a  reducing 
optic  is  required  to  achieve  the  requisite  size.  The  reducing  optic  is  used  to  minify  the 
composite  image  formed  by  the  channel  lens  array  such  that  it  fits  entirely  within 
the  confines  of  the  focal  plane  array.  The  reducing  optic  must  be  large  enough 
(compared  to  the  size  of  the  focal  plane  array)  to  ensure  that  the  reducing  optic 
does  not  inadvertently  define  the  system  field  stop. 

The  field  stop  defines  the  limits  of  the  target  scene  for  each  channel  image 
on  the  FPA.  As  shown  in  figure  D.3,  the  channel  stop  is  implemented  in  the  target 
plane.  This  stop  can  be  realized  physically  by  either  placing  an  aperture  in  front 
of  the  illuminating  source  or  by  placing  an  aperture  between  the  target  and  lens 
array  (though  still  in  the  target  plane).  The  former  option  allows  for  either  front 
or  back  illumination  of  the  target  and  is  therefore  the  preferred  method.  Target 
illumination  is  provided  by  an  LED  source.  Though  LEDs  are  highly  directional, 
target  illumination  will  be  uniform  so  long  as  the  solid  angle  defined  by  the  limiting 
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aperture  is  small  compared  to  the  LED  viewing  angle  (which  is,  by  convention,  the 
angle  between  half  intensity  points  on  either  side  of  direct  viewing). 

D.2.1  System  Design  Restrictions.  Despite  the  relatively  small  number 
of  components  involved,  there  are  a  number  of  important  restrictions  considered 
during  the  design  process.  These  restrictions  include  the  location  and  size  of  the 
system  stops,  the  size  and  spacing  of  the  composite  image  on  the  FPA,  and  the 
space  available  in  the  lab  for  layout  of  the  system  end  to  end.  The  following  section 
describes  each  of  these  restrictions  in  detail. 

As  mention  previously,  each  stop  in  the  channel  array  is  intended  to  serve  as 
an  aperture  stop.  The  image  of  the  optical  element  in  object  space  that  subtends 
the  smallest  angle  at  the  on-axis  position  of  the  target  is  the  entrance  pupil.  The 
physical  element  that  corresponds  to  the  entrance  pupil  is  the  aperture  stop  [31]. 
This  criterion  is  satisfied  when: 


Ds  < 


fc 

d-fc 


Dro 


o  — 


dfc 

d-fc 


(D.6) 


where  Dr  and  Ds  are  the  diameters  of  the  reducing  lens  and  channel  stops,  fc 
is  the  focal  length  of  the  channel  lenses,  and  o  and  d  are  as  defined  in  figure  D.3. 
In  actual  application  this  criterion  is  almost  always  met  because  the  stops  in  the 
channel  array  are  of  a  much  smaller  diameter  than  the  reducing  lens.  The  upper 
limit  on  Ds  is  the  width  of  a  channel  lens. 


Additionally,  the  focal  plane  array  must  serve  as  the  composite  image  field  stop. 
The  image  in  object  space  of  the  optical  element  that  subtends  the  smallest  angle 
at  the  on  axis  position  of  the  entrance  pupil  is  the  entrance  window.  The  physical 
element  that  corresponds  to  the  entrance  window  is  the  field  stop  [31].  Consequently, 
the  reducing  lens  meet  the  following  condition: 


113 


(D.7) 


n  d  WFPA 

r>o^r 

where  WFpa  is  the  width  of  the  focal  plane  array,  M  is  the  overall  system 
magnification,  and  Dr  is  the  diameter  of  the  reducing  lens.  Implicit  in  equation  (D.7) 
is  the  assumption  that  the  channel  array  is,  in  fact,  the  location  of  the  aperture  stop. 

The  sensor  design  requires  that  all  channels  be  imaged  simultaneously  on  a 
single  focal  plane  array.  As  shown  in  figure  D.5,  several  new  quantities  are  required 
to  mathematically  define  this  requirement  in  addition  to  those  already  defined:  the 
width  of  the  field  stop,  Dfs,  and  the  magnification  of  the  channel  lens  alone,  Mc. 


Figure  D.5:  The  projected  image  on  the  FPA. 

Figure  D.5  demonstrates  two  obvious  restrictions  on  the  optical  setup.  First 
the  longest  dimension  of  the  composite  image  must  be  less  than  the  width  of  the 
focal  plane: 


MM~lDc  +  MDfs  <  WFPA  (D.8) 

where  Dc  is  the  diameter  of  a  channel  lens  (which  is,  in  general,  different  than  the 
diameter  of  the  aperture  stop).  Second,  the  images  themselves  must  not  overlap: 
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^  >  Mc  (D.9) 

Ufs 

Ideally,  Mc  should  be  made  as  large  as  possible  to  ensure  a  maximum  number  of 
samples  across  the  target  image  (i.e.  the  most  efficient  use  of  available  detectors  on 
the  focal  plane  array.) 

The  sensor  as  described  is  primarily  intended  to  be  a  lab  instrument  therefore 
compact  packaging  requirements  are  not  as  restrictive  as  those  for  held  instruments. 
Regardless,  lengthy  optical  systems  are  more  difficult  to  align  and  lab  space  is  often 
lacking.  As  such,  overall  system  length,  L ,  is  an  important  design  criterion. 

L  =  o  +  d  +  i  (D.10) 


where  o,  d,  and  i  are  defined  in  figure  D.3.  Given  the  requirement  for  a  specific 
system  magnification  and  an  object  distance,  both  i  and  d  can  be  solved  for  by 
repeated  application  of  the  thin  lens  equation  as  shown  in  [16]: 


M  = 
i  = 


fci 


d{o—fc)—ofc 

fcfrO 


frd- 


O-fc 


d-fr—^r 
1  O-fc 


(D.ll) 


where  the  only  thing  yet  to  be  defined  is  fr,  the  focal  length  of  the  reducing  lens. 

Within  these  confines,  it  is  best  to  maximize  the  total  system  magnification 
(once  the  required  held  stop  diameter  is  determined)  to  ensure  a  maximum  number 
of  pixels  across  the  target.  This  step  should  not  be  confused  with  the  PSF  sampling 
requirements  from  D.5,  which  are  still  necessities.  Recall  that  PSF  sampling  depends 
on  the  radius  of  the  aperture  stop  and  can  be  made  quite  small  at  the  cost  of  reduced 
throughput.  Table  D.2  contains  the  design  used  to  produce  the  laboratory  example 
in  Chapter  IV. 
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Table  D.2:  A  workable  parallel  polarimetric  imager 


fc 

Dc 

Ds 

fr 

Dr 

25  cm 

12.5  mm 

2.38  mm 

10  cm 

5.08  cm 

0 

i 

d 

Dfs 

Ax 

55  cm 

8  cm 

5  cm 

2.5  cm 

16  um 

D.2.2  Advantages  and  Disadvantages.  The  primary  advantage  of  this  sys¬ 
tem,  which  should  not  be  understated,  is  the  fact  that  all  channels  image  simul¬ 
taneously.  The  primary  disadvantages  are  added  complexity  and  reduced  field  of 
view  (when  compared  to  a  single  channel  imager)  and  the  necessity  of  placing  the 
field  stop  at  the  target.  This  design  could  readily  be  improved  by  incorporating  a 
field  stop  into  the  optical  system.  To  do  so,  a  stopped,  intermediate  image  could  be 
formed  at  the  location  of  the  target  in  the  current  design. 
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